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" /  Recent  experiments  have  shown  that  when  the  acoustic  driving 
frequency  is  near  one  of  ^-Ehe  bubble's  harmonic  resonances,  the 
theoretical  values  predicted  by  the  Rayleigh-Plesset  equation  are 
inconsistent  with  observed  values.  This  inconsistency  lead  Prosperetti 
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different  acoustic  cavitation  theories. 


>  XI  c 


A'-;. 


r\> 


•  rs 


Accession  For 

NT IS  GRA&I 
DTIC  T.'R 
Unannounced 
Jjclif  l  c.'it  i  on.. 


S/N  C  102-  LF-  on-  6601 


Unclassified 


SECURITY  CLASSIFICATION  OF  THIS  PACEfHTi.n  O.f.  Cnt.,,d) 


TABLE  OF  CONTENTS 


Page 


LIST  OF  FIGURES .  vi 

INTRODUCTION  .  1 

Chapter 

1.  THEORETICAL  DEVELOPMENT  .  6 

A.  The  acoustic  cavitation  equations  .  6 

B.  An  exact  formulation  for  the  internal  pressure  ....  15 

C.  Limitations  of  the  exact  formulation  .  21 

2.  NUMERICAL  SOLUTION  OF  THE  EQUATIONS  .  25 

A.  Dimensionless  form  of  the  equations .  25 

B.  Finite  difference  equations  .  30 

C.  Error  analysis  .  . .  38 

D.  Applicability .  42 

3.  COMPUTATIONAL  RESULTS .  44 

A.  Radius  verses  time  curves . 44 

B.  Frequency  response  curves  ....  .  52 

C.  Levitation  numbers  .  ........  61 

D.  Thermodynamics  of  the  bubble  interior .  67 

4.  A  SYSTEMATIC  APPROACH  TO  CHAOTIC  BUBBLE  MOTION  .  81 

5.  CONCLUSIONS  AND  TOPICS  FOR  FUTURE  STUDY  .  87 

REFERENCES .  92 


V 


Chapter  Page 

APPENDIX .  97 

A.  Computer  program  for  the  exact  formulation  .  97 

R.  Computer  program  for  the  Ray leigh-Plesset 

Polytropic  equation  ....  .  106 

BIOGRAPHICAL  SKETCH  OF  THE  AUTHOR  .  112 


LIST  OF  FIGURES 


Figure  Page 

1.  Thermal  conductivity  of  air  as  a 

function  of  temperature  .  29 

2.  Flowchart  for  acoustic  cavitation  program  .  37 

3.  Normalized  bubble  radius  as  a  function  of  time 
for  the  analytic  first  order  solution  and  our 

new  exact  formulation  .  40 

4.  Normalized  bubble  radius  as  a  function  of  time 
for  the  Rayleigh-Plesset  Polytropic  equation 

and  the  exact  formulation .  41 

5.  Normalized  bubble  radius  as  a  function  of  time  for 

the  RPP,  RPE,  and  CE  methods.  R  -50  microns,  Pa- 

0.4  bar,  and  f-0.4f .  46 

’  o 

6.  Normalized  bubble  radius  as  a  function  of  time  for 

the  RPP  and  CE  methods.  R  -50  microns,  Pa-0.4  bar, 

and  f-0.8f  ° .  47 

o 

7.  Normalized  bubble  radius  as  a  function  of  time  for 

the  RPP  and  CE  methods.  R  *50  microns,  Pa-0.6  bar, 

and  f-0.8f  0 .  48 

o 

8.  Normalized  bubble  radius  as  a  function  of  time  for 

the  RPP  and  CE  methods.  R  -10  microns,  Pa-0.6  bar, 

and  f -0.4 f  . ° .  50 

o 

9.  Normalized  bubble  radius  as  a  function  of  time  for 

the  RPP  and  CE  methods.  R  -10  microns,  Pa-0.6  bar, 

and  f-0.8f  ° .  51 

o 

10.  Maximum  normalized  bubble  radius  as  a  function  of 
normalized  frequency  for  the  CE  method.  R  -10 

microns.  Pa-  0.3,  0.5,  and  0.7  bar  .  .  .  .  53 

11.  Maximum  normalized  bubble  radius  as  a  function  of 
normalized  frequency  for  the  CE  method.  R  -50 

microns,  Pa-  0.2,  0.3,  0.4,  and  0.5  bar  .° .  54 


BwKH 


»  '  -  '/  V  /  • 

’  .  *  VT  •  *  «  V 


Vi 


vii 


Figure  page 

12.  Maximum  normalized  bubble  radius  as  a  function  of 
normalized  frequency  for  the  CE  method.  R  *100 

microns,  Pa*0.1,  0.2,  0.3,  0.4,  0.5,  0.6,  and  0.7  bar  .  .  55 

13.  Maximum  normalized  bubble  radius  as  a  function  of 
normalized  frequency  for  the  CE  method.  R  =10 
microns,  Pa=0.7  bar,  asterisks  denote  points  from 

reference  [52] .  57 

14.  Maximum  normalized  bubble  radius  as  a  function  of 
normalized  frequency  for  the  CE  method.  R  =100 
microns,  Pa»0.7  bar,  asterisks  denote  points  from 

reference  [17] .  58 

15.  Maximum  normalized  bubble  radius  as  a  function  of 
normalized  frequency  for  the  CE  method.  R  =10 
microns,  Pa*0.7  bar,  asterisks  denote  ultraharmonic 


resonance  peaks  from  reference  [52] . .  .  60 

16.  Levitation  number  as  a  function  of  the  normalized 
radius  for  the  CE  and  RPP  methods.  Freq .=22 .2kHz , 

Pa=0.095  bar,  asterisks  denote  data  from  ref.  [16]  ....  64 

17.  Levitation  number  as  a  function  of  the  normalized 
radius  for  the  CE  and  RPP  methods.  Freq .=22 ,2kHz , 

Pa=0.155  bar,  asterisks  denote  data  from  ref.  [16]  ....  65 


18.  Levitation  number  as  a  function  of  the  normalized 
radius  for  the  CE  and  RPP  methods.  Freq .=22. 2kHz, 

Pa*0.190  bar,  asterisks  denote  data  from  ref.  [16]  ....  66 

19.  Temperature  profile  for  a  1  micron  bubble  driven 
at  2.1  bar  and  0.4  times  the  linear  resonance 
frequency.  Values  for  the  minimum  and  maximum 

radius  are  shown .  69 

20.  Temperature  profile  for  a  10  micron  bubble  driven 
at  0.9  bar  and  0.4  times  the  linear  resonance 
frequency.  Values  for  the  minimum  and  maximum 

radius  are  shown .  70 

21.  Temperature  profile  for  a  100  micron  bubble  driven 
at  0.8  bar  and  0.4  times  the  linear  resonance 
frequency.  Values  for  the  minimum  and  maximum 

radius  are  shown . 71 


^Tv 


viii 


Figure 

22.  Center  temperature  as  a  function  of  time  for  a  1 
micron  bubble  driven  at  1.4  bar,  and  frequencies 

of  0.4  and  0.8  times  the  linear  resonance  frequency  .  .  . 

23.  Center  temperature  as  a  function  of  time  for  a  10 
micron  bubble  driven  at  0.6  bar,  and  frequencies 

of  0.4  and  0.8  times  the  linear  resonance  frequency  .  .  . 

24.  Center  temperature  as  a  function  of  time  for  a  100 
micron  bubble  driven  at  0.45  bar,  and  frequencies 

of  0.4  and  0.8  times  the  linear  resonance  frequency  .  .  . 

25.  Internal  pressure  as  a  function  of  time  for  a  1 
micron  bubble  driven  at  1.4  bar,  and  frequencies 

of  0.4  and  0.8  times  the  linear  resonance  frequency  .  .  . 

26.  Internal  pressure  as  a  function  of  time  for  a  10 
micron  bubble  driven  at  0.6  bar,  and  frequencies 

of  0.4  and  0.8  times  the  linear  resonance  frequency  .  .  . 

27.  Internal  pressure  as  a  function  of  time  for  a  100 
micron  bubble  dirven  at  0.45  bar,  and  frequencies 

of  0.4  and  0.8  times  the  linear  resonance  frequency  .  .  . 

28.  Radius,  internal  pressure,  and  center  temperature 

as  a  function  of  time  for  a  50  micron  bubble  driven 
at  0.6  bar  and  a  frequency  of  0.4  times  the  linear 
resonance  frequency  ........  .  ... 


Page 


73 


74 


75 


77 


78 


79 


80 


29.  Feigenbaum  tree  for  the  RPP  equation.  Rq=50 

microns  and  f/f  =0.5..  .......... 

o 

30.  Feigenbaum  tree  for  the  CE  equations.  Rq=50 

microns  and  f/f  *0.5 . . 

o 


85 

86 


w.-h.  . 


INTRODUCTION 


The  focus  of  this  dissertation  is  on  numerical  modeling  of  stable 
acoustic  cavitation.  Neppiras  [1]  defines  stable  acoustic  cavitation  as 
the  oscillation  of  cavities  or  bubbles  about  some  equilibrium  size.  In 
general  the  oscillations  need  not  be  linear,  but  may  vary  with  time  in  a 
very  complex  manner.  The  adjective  stable  implies  that  the  bubble 
exists  for  many  cycles  without  breaking  up  or  dissolving.  On  the  other 
hand,  a  transient  cavitation  event  is  one  which  usually  exists  for  one 
cycle  or  less.  The  bubble  grows  to  several  times  its  equilibrium  size 
and  then  collapses  violently  breaking  into  many  smaller  bubbles.  We 
will  only  deal  with  stable  cavitation  here. 

Before  pressing  onward  with  the  latest  acoustic  cavitation 
theories,  let  us  briefly  examine  the  short  history  of  this  field  of 
Physics.  Interest  in  cavitation  dates  back  to  Besant  [2]  in  1859. 
However,  there  was  no  significant  theoretical  work  in  this  field  until 
that  of  Rayleigh  [3]  in  1917.  This  paper  by  Rayleigh  described  the 
collapse  of  a  spherical  cavity.  Of  course  Lord  Rayleigh  is  most  famous 
for  his  treatise  The  Theory  of  Sound  which  he  started  writing  while  on 
vacation  in  Egypt  in  1872.  This  monumental  work  took  five  years  to 
reach  the  press  and  it  is  a  credit  to  Rayleigh's  abilities  that  it  is 
still  considered  today  as  the  foundation  of  physical  acoustics. 

Between  Rayleigh's  paper  in  1917  and  the  late  1940's,  all 
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theoretical  work  in  cavitation  was  concerned  with  hydrodynamically- 
generated  cavitation,  usually  by  propellers  on  ships.  In  1949  Blake  [4] 
published  the  first  systematic  study  on  acoustic  cavitation.  Shortly 
after  this  publication,  the  classic  papers  by  Plesset  [5],  Neppiras  and 
Noltingk  [6,7],  and  Poritsky  [8]  appeared  that  significantly  advanced 
the  theory.  It  was  about  this  time  that  many  groups  began  working  on 
the  problem  of  acoustic  cavitation.  Plesset  added  variable  pressure  and 
surface  tension  terms  to  the  theory  of  Rayleigh  and  for  his  contribution 
the  resulting  equation  with  a  damping  term  added  by  Poritsky  is  called 
the  Rayleigh-Plesset  equation.  This  equation  is  derived  and  discussed 
in  chapter  1. 

The  basic  problem  of  acoustic  cavitation  is  to  find  the  pressure 
and  velocity  fields  of  the  two  phase  medium  consisting  of  a  bubble  or 
cavity  and  a  surrounding  fluid  which  theoretically  extends  to  infinity. 
Recent  advances  in  ultrasonic  instruments  used  for  medical  purposes  have 
made  a  knowledge  of  the  internal  temperature  of  an  oscillating  bubble 
important  as  well.  High  temperatures  in  the  interior  of  the  bubble  can 
produce  free  radicals  which  could  be  dangerous  to  biological  systems 
[9].  In  practice  the  bubble's  radius  as  a  function  of  time  is  also  of 
interest  as  well.  Although  the  new  theoretical  model  is  very  complex 
and  requires  a  lot  of  computer  time  for  its  solution,  all  of  the 
quantities  mentioned  above  can  be  obtained  from  its  solution. 

In  this  study  numerical  methods  are  used  extensively  to  solve  the 
basic  equations  of  acoustic  cavitation  for  the  case  of  a  spherically 
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oscillating  bubble  in  which  the  internal  pressure  is  defined  by  a  method 
suggested  by  Prosperetti  [10].  A  thorough  comparison  with  previous 
results  will  be  made  to  see  just  what  parameter  values  necessitate  the 
use  of  this  more  complex  model.  In  chapter  4  it  will  be  shown  that  this 
new  theory  leads  to  a  different  type  of  coupling  between  the  liquid  and 
gas  phases  through  the  use  of  "Feigenbaum  trees"  [11,12]. 

The  most  elementary  model  of  acoustic  cavitation  considered  here  is 
the  Rayleigh-Plesset  equation,  derived  by  Plesset  [5]  with  the  viscous 
damping  term  added  by  Poritsky  [8].  Extensive  numerical  investigations 
were  carried  out  for  this  equation  by  Noltingk  and  Neppiras  [6,7]  and 
later  by  Flynn  [13,14].  Since  that  time  several  approximate  analytical 
solutions  to  the  Rayleigh-Plesset  equation  have  been  derived,  including 
a  solution  to  second  order  in  the  asymptotic  expansion  by  Prosperetti 
[15].  These  solutions  were  thought  to  be  adequate  until  recent 
experimental  data  showed  a  large  deviation  between  theory  and  experiment 
for  driving  frequencies  in  the  vicinity  of  a  harmonic  of  the  resonance 
frequency  of  the  bubble.  It  will  be  shown  that  the  new  theory  is  closer 
to  experimental  data  In  these  sensitive  frequency  ranges  and  thus  is  a 
better  model  for  the  thermal  damping  in  the  bubble. 

The  first  chapter  contains  a  brief  derivation  of  some  of  the  basic 
equations  in  cavitation.  Included  in  this  section  are  the  well-known 
Rayleigh  equation  and  the  Rayleigh-Plesset  equation.  The  more  advanced 
equations  for  radial  bubble  oscillations  are  discussed  as  well  and 
relationships  between  them  are  pointed  out.  After  establishing  the 
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heating  term  is  given  by  Prosperetti  [8]  as 
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where  T  is  the  internal  temperature,  the  specific  heat  of  the  gas  at 
constant  pressure,  and  K  is  the  thermal  conductivity  of  the  gas. 
Multiplication  of  equation  (38)  by  c^T  and  adding  it  to  equation  (40) 
results  in 
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Using  the  ideal  gas  assumption  it  is  easy  to  show  that 
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Sustituting  equations  (42)  and  (37)  into  (41)  results  in 
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Because  of  the  previously  assumed  spherical  symmetry,  one  can  multiply 
both  sides  of  equation  (43)  by  dV  and  integrate.  Applying  the 
divergence  theorem  to  the  second  term  yields  the  following 
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Carrying  out  the  integration  yields 
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Now,  solving  for  the  gas  velocity  U  gives 


where  p,  P,  U,  and  t  denote  the  density,  pressure,  radial  velocity,  and 
viscous  stress  tensor  of  the  gas  in  the  bubble.  r  =  0  corresponds  to 
the  center  of  the  bubble  and  r  =  R  corresponds  to  the  wall  of  the  bubble 
of  radius  R.  The  boundary  conditions  for  the  velocity  are 

U(r=0,t)  =  0  and  U(r=R,t)  =  dR/dt.  (36) 

From  the  assumption  of  spherically  symmetric  oscillations  and 
radial  gas  velocities,  it  is  clear  that  the  internal  pressure  is  only 
a  function  of  the  radial  coordinate  r  and  time  t.  In  section  C  of  this 
chapter  we  will  show  that  order  of  magnitude  estimates  reveal  that  we 
can  replace  equation  (35)  by 

P1  =  PA(t)  or  P£  =  Pi(R)  (37) 

since  R  *  R(t). 

Next,  the  continuity  equation  is  considered  (  see  for  example 
Batchler  [23]  ) 
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Finally,  the  conservation  of  energy  equation  without  the  viscous 
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B.  An  exact  formulation  for  the  internal  pressure 

The  need  for  a  more  accurate  expression  for  the  internal  pressure 
of  a  bubble  as  a  function  of  both  the  bubble  radius  and  time  led 
Ijosperetti  [10]  to  a  direct  approach  using  the  conservation  equations 
for  the  interior  of  a  bubble.  The  derivations  that  follow  are  along  the 
lines  given  by  Prosperetti. 

Certain  assumptions  about  the  geometry  of  the  bubble  oscillations 
and  the  nature  of  the  gases  and  liquids  will  be  made  from  the  start  of 
the  derivation.  These  assumptions  set  some  limitations  on  the  theory 
that  follows.  These  limitations  will  be  discussed  in  detail  in  section 
C  of  this  chapter. 

The  first  assumption  made  is  that  the  bubble  oscillations  are 
spherical  and  that  the  velocity  of  the  gas  in  the  bubble  is  only  a 
function  of  a  radial  coordinate  and  time.  Fanelli  et.  al.  [41]  have 
shown  that  diffusion  of  gas  in  and  out  of  the  bubble  is  only  important 
at  low  ambient  pressures  and  hence  the  assumption  of  no  mass  transport 
across  the  bubble  wall  is  appropriate.  This  study  is  only  concerned 
with  bubble  oscillations  in  water  at  or  near  room  temperature.  Since 
the  vapor  pressure  of  water  at  room  temperature  is  small  compared  to 
that  of  the  gas  in  the  bubble,  it  is  neglected  in  our  derivation.  The 
final  assumption  is  that  the  gas  in  the  bubble  is  an  ideal  gas. 

The  momentum  conservation  equation  can  be  written  as  follows  (  see 


for  example  Aris  [24]) 
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Until  recently  the  polytropic  approximation  was  thought  to  be  quite 
accurate  over  a  wide  range  of  driving  pressures,  frequencies,  and  bubble 
radii  that  are  of  interest.  However,  Crum  and  Prosperetti  [19]  have 
shown  that  the  polytropic  approximation  can  give  quite  different  results 
when  the  driving  frequency  is  in  the  vicinity  of  one  of  the  harmonics  of 
the  resonance  frequency  of  the  bubble.  It  is  this  difficulty  that  led 
Prosperetti  to  search  for  a  better  expression  for  the  internal  pressure, 
which  is  the  subject  of  the  next  section. 
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reference  [18].  Even  though  equation  (28)  looks  quite  different  from 
(26),  substitution  of 

K  =  (P  -  ?J/Pm  (32) 

into  (28)  with  C  held  constant  simplifies  (28)  to  (26).  Thus,  these  two 
equations  are  actually  quite  similar. 

To  this  point,  only  the  fundamental  equations  for  radial 
oscillations  have  been  discussed  and  not  the  problem  of  how  to  model  the 

internal  pressure.  For  some  time  now  the  internal  pressure  has  been 

approximated  by  an  expression  involving  a  polytropic  exponent  of  the 
form  given  by  Zwick  [40] 

P±  =  Po(Ro/R)3K  .  (33) 

Here  k  is  the  polytropic  exponent  and  Pq  is  the  internal  pressure  of  the 
bubble  at  equilibrium  defined  by 

Po  =  P^  +  2a/R0  ,  (34) 

where  Rq  is  the  equilibrium  radius  of  the  bubble  (  note  that  equation 
(30)  is  of  this  form  as  well  ).  Using  this  approximation  for  P^.,  it  is 
straightforward  to  solve  equation  (25)  numerically  using  a  fourth  order 
Runge-Kutta  method  to  a  high  degree  of  precision.  There  are  of  course 
other  numerical  methods  to  solve  this  type  of  second  order  nonlinear 
ordinary  differential  equation,  but  the  Runge-Kutta  method  is  probably 
the  most  widely  used  and  suffices  here. 
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by  the  condition  on  the  normal  stresses 


Pt(R,t)  «  P(R,t)  +  2o/R  +  4pR/R.  (27) 

It  is  easy  to  show  that  as  c  approaches  infinity,  equation  (26) 
simplifies  to  the  Rayleigh-Plesset  equation  (25). 

There  are  numerous  equations  equivalent  to  or  similar  to  equation 
(26)  in  the  literature  [30,32,35,37],  Prosperetti  [38]  has  discussed 
this  fact  and  argues  in  favor  of  equation  (26)  for  reasons  of  better 
numerical  convergence.  For  this  reason  equation  (26)  is  used  as  one  of 
the  fundamental  equations  for  this  study. 

Another  equation  which  is  very  popular  in  the  literature  is  that  by 
Gilmore  [30];  however,  the  form  used  here  is  given  by  Lauterborn  and 
Suchla  [18]  as 


(28) 


where  H  is  the  free  enthalpy,  which  for  water  is  given  by 

H  ,  AL-  { [P(R)  +  B]  (n-1)y,n  _  (p  +  B)  (29) 

n  -  1  po 

The  pressure  at  the  bubble  wall  P(R)  and  the  speed  of  sound  c  are  given 
by 


P(R) 


[».  ♦  f 1 

V 

Rj 

rJ 

3Y 


2a  4pR 
R  “  R 


(30) 


and 


C  «  [C§  +  (n-l)H]  . 


(31) 


Typical  values  of  the  constants  in  these  equations  are  given  in 
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time,  the  solution  of  (21)  is  of  the  form 

k  -k 
♦  =  Ar  +  Br  . 


(23) 


2 

The  boundary  conditions  imply  that  A  -  0,  k  =  1,  and  B  =  -RR  .  Thus, 
equation  (23)  becomes 


$ 


(24) 


Now  substituting  equations  (24)  and  (20)  into  (22)  yields  the  Rayleigh- 
Plesset  equation 


RR  +  |R 


yR,t)  -  ps(t>  -  ¥-■¥ 


(25) 


There  are  many  equations  in  the  literature  describing  spherical 
bubble  oscillations  [37].  Probably  the  best-known  of  these  equations  is 
the  Rayleigh-Plesset  equation  (25).  One  of  the  limitations  of  this 
equation  is  the  assumption  of  an  incompressible  liquid.  This  implies 
that  the  speed  of  sound  in  the  liquid  must  be  infinite.  Prosperetti 

[38]  has  shown  that  following  the  procedure  used  by  Keller  and  Miksis 

[39]  the  following  equation  of  motion  for  the  radial  bubble  oscillations 
can  be  obtained 


•  <4 

1  -  * 

3*2 
RR  +  4R2 

f  •  \ 

1  -  *  = 

f  • 

1  +  * 

CJ 

2 

3c 

c 

P(R,t)  -  Ps(t+R/c) 
P 


R  dP(R,t) 
pc  dt 


(26) 


Here  P(r,t)  is  the  pressure  on  the  external  side  of  the  bubble  wall. 
The  pressure  P(R,t)  is  related  to  the  internal  bubble  pressure  P  (r,t) 
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!£  +  {(V*)2  +  h  =  0. 


(18) 


Setting  the  right  hand  side  of  equation  (18)  equal  to  zero  implies  that 
p  approaches  zero  as  r  approaches  infinity.  The  other  boundary 
conditions  on  equations  (17)  and  (18)  are 


9  (ft 

& 


r-R(t) 


(19) 


and 

P  (R, t)  -  P(R.t)  +  4uR/R,  (20) 

8  K 

where  P  ia  the  pressure  in  the  bubble,  P  the  pressure  on  the  liquid 

8 

side  of  the  interface,  and  o  the  surface  tension  of  the  liquid. 
Equations  (17)  and  (18)  together  with  the  boundary  conditions  (19)  and 
(20)  have  been  the  subject  of  considerable  theoretical  work  [30-37]. 

The  equations  above  can  be  simplified  by  assuming  in  the  continuity 
equation  that  the  speed  of  sound  c  is  infinite,  implying  that  the 
density  0  ia  constant.  These  assumptions  reduce  equation  (17)  to 
Laplace's  equation 


7*4  -  0 


(21) 


and  equation  (18)  to 


dp 

at 


+ 


+ 


o 


P«(t) 


(22) 


where  P  is  the  variable  part  of  the  pressure  in  the  liquid.  Since 
s 

these  equations  are  only  dependent  on  the  radial  space  coordinate  and 
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t 


where 


(13) 


Here,  p  is  the  viscosity  of  the  liquid  and  p  is  the  density  of  the 
liquid,  both  assumed  to  be  constant.  For  our  purposes  the  term  in 
equation  (12)  involving  the  viscosity  is  small,  and  together  with 
equation  (9)  we  can  write  equation  (12)  as 


3u  3u  _  J_3P 

3t  3r  =  p3r  ' 


(14) 


Considered  here  are  bubble  oscillations  where  the  internal 
temperature  remains  low  enough  so  that  evaporation  and  condensation  at 
the  interface  can  be  neglected.  This  assumption  allows  us  to  treat  the 
liquid  as  isothermal  and  an  equation  of  state  involving  the  pressure  and 
density  can  be  assumed,  which  leads  to  the  following  expressions 
involving  the  molar  enthalpy  h  and  the  speed  of  sound  c  [28] 

pdh  *  dP  and  c2dp  =  dP.  (15) 


Assuming  a  velocity  potential  <£  such  that 


U  =  V<f)  and  u  *  j^-, 
3r 


(16) 


the  continuity  equation  (11)  becomes 


v2$  +  -2 

c 


3h  5h 

3t  3r  3r 


=  0. 


(17) 


Performing  one  integration  results  in  the  Bernoulli  integral 


V-'.Cis.v.'-s.V 
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This  famous  equation  can  be  very  useful  as  a  starting  point  in  deriving 
acoustic  cavitation  theories. 

Next,  the  acoustic  cavitation  equations  are  derived  by  starting 
from  first  principles  -  namely  the  conservation  equations.  Since 
spherically  symmetric  oscillations  are  dealt  with  exclusively,  the 
velocity  in  spherical  coordinates  is  given  by 

U  *  u^r*  where  u  *  |u|.  (9) 


There  are  numerous  texts  in  fluid  mechanics  and  transport  phenomena,  but 
several  advanced  texts  that  are  particularly  useful  are  those  by 
Batchelor  [23],  Aris  [24],  Bird,  Stewart,  and  Lightfoot  [25],  and  Landau 
and  Lifshitz  [26].  Some  intermediate  texts  that  are  useful  are  those  by 
Lu  [27],  John  [28],  and  Bertin  and  Smith  [29J. 

The  continuity  equation  or  conservation  of  mass  equation  can  be 
written  as 


H  +  v  .  (pU) 


Using  equation  (9)  for  U  this  becomes 


3P  +  U3_P  +  £  L_ 
3t  uar  r23r 


r2u 


0  , 


(10) 


(ID 


where  r  is  the  radial  coordinate  measured  from  the  center  of  the  bubble. 
The  equation  of  motion  or  conservation  of  linear  momentum  equation  known 
as  the  Navier-Stokes  equation  can  be  written  as 


UU  7 

PfT  58  -VP  +  ViV2U  , 


(12) 


k  ■  *■  <kili 


•  *  >  '  «  *  •  *  .  *  .  •  .  ,  •  . 


-  ■£-  if- 
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u  =  R2R/r2 


(3) 


Substituting  equation  (3)  into  (2)  and  integrating  results  in  the 
kinetic  energy  given  by 


K  -  2 7T p R  3  R2  . 


(4) 


The  power  J  can  be  defined  by 


(5) 


where  F  is  the  force  which  equals  pressure  times  area  and  ii  is  the 
velocity  of  the  fluid.  Since  the  force  and  velocity  are  both  radially 
directed  and  antiparallel,  the  power  at  infinity  becomes 


J  *  1  im 
r->°° 


f-p  | 

.  J  l  4 


(6) 


Substituting  equation  (3)  for  u  and  recalling  that  the  pressure  at 


infinity  is  PQ  we  get 


J  ■  -P  4ttR2R. 
o 


The  time  derivative  of  the  kinetic  energy  (4)  is  given  by 


J  =  2irp  (3RZR3  +  2R3RR) 


(7) 


(8) 


Equating  the  power  (7)  with  the  time  derivative  of  the  kinetic  energy 
(8)  results  in  the  Rayleigh  equation 


RR  +  |r2  =  -P  /p. 


(1) 
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Chapter  1 

Theoretical  development 

A.  The  acoustic  cavitation  equations 

Before  discussing  the  derivation  of  the  acoustic  cavitation 
equations,  let  us  first  recall  the  method  of  obtaining  the  Rayleigh 
equation  [3] 

RR  +  |r2  -  -Po/p.  (1) 

Here,  R,  PQ,  and  p  denote  the  bubble  radius,  pressure  at  infinity,  and 
density  of  the  liquid  respectively.  Dots  denote  differentiation  with 
respect  to  time.  It  is  assumed  that  the  liquid  is  inviscid, 
incompressible,  and  has  no  surface  tension;  also  the  internal  pressure 
is  assumed  to  be  zero. 

Following  the  method  of  Batchelor  [23],  the  time  rate  of  change  of 
kinetic  energy  of  the  fluid  is  assumed  to  be  equal  to  the  power 
developed  at  infinity.  The  kinetic  energy  of  the  fluid  is  defined  by 

K  *  Hp  4irr  2u  zdr ,  (2) 

Jr 

where  u  is  the  magnitude  of  the  fluid  velocity.  The  continuity  equation 
together  with  the  above  assumptions  imply  that  the  magnitude  of  the 
velocity  can  be  written  as 
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The  temperature  profixts  will  prove  very  helpful  in  determining  what 
type  of  damage  could  be  caused  by  therapeutic  and  diagnostic  ultrasound 
instruments . 

The  results  in  chapter  4  are  separated  from  those  in  chapter  3 
because  they  come  from  a  completely  different  approach  to  understanding 
how  the  model  behaves  under  certain  initial  conditions.  The  system  of 
equations  describing  the  bubble  motion  is  viewed  as  a  dynamical  system. 
In  analyzing  the  dynamical  system  some  of  the  recent  techniques  such  as 
those  by  Ott  [20],  Huberman  and  Crutchfield  [21],  and  Grebogi,  Ott,  and 
Yorke  [22]  are  employed.  This  analysis  will  show  that  the  type  of 
damping  used  in  the  equations  dramatically  alters  the  results  of  these 
techniques. 

Finally,  in  chapter  5  some  general  conclusions  are  drawn  from  the 
information  compiled  and  suggestions  are  made  as  to  what  direction  new 
research  should  follow.  Ideas  for  both  new  experimental  work,  not 
discussed  in  previous  chapters,  and  improvements  in  the  theoretical 
model  will  be  presented. 


fundamental  equations  of  acoustic  cavitation,  the  internal  pressure  is 
then  considered.  The  method  of  Prosperetti  [10]  will  be  employed  in 
developing  the  theory  used  to  solve  for  the  internal  pressure.  Chapter 
1  concludes  with  a  brief  error  analysis  for  all  the  terms  discarded  in 
the  conservation  equations,  which  are  used  to  solve  for  the  internal 
pressure. 

Chapter  2  deals  with  the  numerical  solution  of  the  system  of 
acoustic  cavitation  equations  that  is  derived  in  chapter  1.  The  first 
step  is  to  find  the  dimensionless  form  of  all  the  equations  that  will  be 
used  in  computer  programs.  Next,  the  numerical  techniques  used  to  solve 
this  coupled  system  of  nonlinear  ordinary  and  partial  differential 
equations  along  with  the  finite  difference  representation  of  the 
equations  are  described  in  detail.  A  brief  error  analysis  for  the 
numerical  methods  used  concludes  this  chapter. 

Chapter  3  contains  results  not  previously  known  and  comparisons 
between  the  new  results  and  those  of  previous  theories.  In  particular, 
radius  verses  time  curves  for  several  models  and  frequency  response 
curves  such  as  those  by  Cramer  [17]  and  Lauterborn  and  Suchla  [18]  are 
compared  to  the  new  theory.  A  good  test  to  determine  how  well  the 
theory  estimates  the  damping  is  to  compare  theoretical  levitation  number 
curves  with  the  experimental  ones  obtained  by  Crum  and  Prosperetti  [19]. 
Chapter  3  concludes  with  a  discussion  of  the  internal  thermodynamics  of 
bubbles.  Both  temperature  profiles  and  pressure  curves  are  presented 
as  a  function  of  radius  or  time  for  a  variety  of  initial  conditions. 


c 


if  3T  1  • 

U  =  y“p  (Y“1)KaT-  3rP 


Using  the  boundary  conditions  (36)  and  evaluating  equation  (46)  at  r  *  R 
gives  the  following  ordinary  differential  equation  for  P 

P  =  |  (y-DKerj  !£  "  yP'R  '  (47) 

'  R  ' 

With  the  approximation  of  uniform  pressure,  equation  (47)  contains 
all  of  the  information  present  in  the  continuity  equation.  Having 
already  simplified  two  of  the  three  original  partial  differential 
equations,  we  turn  our  attention  once  again  to  the  energy  equation. 
Using  the  ideal  gas  assumption,  equation  (40)  becomes 


ySTl  _  p  -  p  .  (KVT), 
Y-ljT(3t  3rj 


where  U  and  dP/dt  are  defined  by  equations  (46)  and  (47)  respectively. 
The  boundary  conditions  for  equation  (48)  are  continuity  of  temperature 
and  heat  fluxes  across  the  bubble  wall.  If  the  liquid  is  cold  enough  so 
that  the  vapor  pressure  of  the  liquid  is  small  compared  to  that  of  the 
gas  (  an  assumption  which  holds  at  room  temperature  ),  one  can  assume 
that  the  temperature  boundary  condition  is 


T(r*R,t)  =  T 


This  assumption  greatly  simplifies  the  solution  of  equation  (48)  and 


will  be  discussed  at  length  in  section  C  of  this  chapter. 

During  the  collapse  phase  of  the  bubble  motion,  the  gas  temperature 
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can  be  exceedingly  high.  For  this  reason  one  must  take  into  account  the 
variation  of  the  thermal  conductivity  K(T)  for  the  gas  in  the  bubble. 
Since  the  temperature  of  the  gas  is  a  function  of  both  position  and 
time,  so  is  the  conductivity.  This  presents  a  problem  on  the  right  hand 
side  of  equation  (48).  The  standard  technique  for  dealing  with  this 
difficulty  is  to  define  a  new  variable  t  as  an  integral  of  the 
conductivity  over  appropriate  limits  such  as 

t  =  [  K(x) dx .  (50) 
JT 

00 

Another  standard  change  of  variables  to  simplify  equation  (48)  is 
to  define  a  fixed  boundary  rather  than  a  varying  one  by  introducing  the 
new  radial  coordinate 

y  -  r  /  R(t).  (51) 
Using  equations  (46),  (50),  and  (51)  along  with  some  algebraic 
manipulations  one  can  write  equation  (48)  as 


3-1  + 

ftJj 

ax 

ax 

' 

at 

[yPR2J 

.3y 

"  y3y 

i. 

|  -  DP  =  , 


(52) 


where  the  Laplacian  is  with  respect  to  y  and  D(P,t)  is  defined  by 


D(P,t) 


K(T) 

Cpp(P,T) 


y-l|K(T)T 
.  Y  j  p 


(53) 


which  is  the  form  of  the  thermal  diffusivity  for  an  ideal  gas.  The 
boundary  conditions  for  t  are  given  by 

t ( y= 1 , t )  =  0.  (54) 

The  substitutions  (50)  and  (51)  are  also  made  in  equations  (46)  and 
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(47).  The  new  forms  of  these  two  equations  add  nothing  to  the  theory 
and  will  be  omitted  here.  They  can  be  found  in  chapter  2,  section  A. 

All  of  the  equations  necessary  to  solve  for  the  bubble's  radius  as 
a  function  of  time  using  the  new  model  for  the  internal  pressure  have 
now  been  derived.  Because  of  the  techniques  employed,  one  also  gets 
P(t)  and  T(r,t)  as  well  as  R(t).  The  temperature  information  will  prove 
to  be  very  important  in  some  practical  problems  discussed  later. 

The  method  of  solution  of  these  equations  is  simple  in  principle. 
The  Rayleigh-Plesset  equation  (25)  or  an  appropriate  equation  containing 
the  compressibility  terms  (26)  is  coupled  to  the  energy  equation  (52), 
the  pressure  equation  (47),  and  U  ■  dR/dt.  One  can  solve  this  system  of 
differential  equations  for  the  radius,  pressure,  velocity,  and 
temperature  as  a  function  of  time  by  using  numerical  methods.  The  exact 
method  of  solution  is  discussed  in  detail  in  chapter  2  section  B. 
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Limitations  of  the  exact  formulation 


In  numerical  modeling  of  physical  problems,  there  is  usually  a 
compromise  between  the  degree  of  accuracy  desired  and  the  amount  of 
computer  power  available  to  do  the  analysis.  As  the  accuracy  of  the 
equations  is  increased  by  adding  smaller  and  smaller  physical  effects, 
it  sometimes  happens  that  a  set  of  equations  is  arrived  at  that  is  too 
complicated  to  solve.  Without  the  aforementioned  approximations  for  the 
internal  pressure  of  a  cavitation  bubble,  the  new  model  would  have  this 
problem.  However,  enough  information  has  been  included  in  the  equations 
to  give  considerable  improvement  over  previous  models  for  the  range  of 
parameter  values  which  are  of  interest.  As  one  ventures  out  of  this 
parameter  range,  one  might  hope,  but  should  not  expect  the  new  theory  to 
agree  very  closely  with  experimental  values. 

The  focus  of  this  study  is  to  show  an  improved  performance  for  the 
fundamental  equations  given  in  section  A  of  this  chapter.  Since  all  of 
these  equations  assume  radial  oscillations,  it  is  assumed  that  this  is 
the  case  for  the  treatment  of  the  internal  pressure.  Prosperetti  [42] 
gives  conditions  favorable  for  non-spherical  oscillations  and  the 
parameter  values  used  exclude  this  type  of  behavior.  Using  a  Levitation 
cell  like  the  one  used  by  Crum  [43],  one  can  easily  detect  the  onset  of 
non-spherical  oscillations.  This  device  is  used  for  experimental 
measurements  and  one  can  determine  if  appreciable  non-spherical 
oscillations  are  occurring  which  might  cause  the  data  to  deviate  from 
those  predicted  by  the  theory.  One  can  conclude  that  any  differences 
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between  theoretical  and  experimental  values  must  be  due  to  some  other 
approximations  we  have  made. 

In  the  derivation  of  the  fundamental  equations  in  section  A  of  this 
chapter  it  was  assumed  that  the  bubble  wall  velocity  is  equal  to  the 
liquid  velocity  at  r  «  R.  Prosperetti  [44]  has  shown  that  this  is  valid 
only  if  little  or  no  mass  transfer  takes  place  across  the  interface. 
The  error  due  to  this  effect  is  of  the  order  of  the  ratio  of  the  density 
of  the  gas  in  the  bubble  to  the  density  of  the  surrounding  liquid.  A 
typical  value  for  this  ratio  is  0.001  or  less. 

Previously,  the  momentum  equation  (35)  was  shown  to  imply  that 
P^  =  P  (t).  This  required  the  assumption  that  all  of  the  terms  in  the 
momentum  equation  contribute  to  but  a  small  spatial  pressure  difference 
in  the  bubble.  Prosperetti  [10]  gives  a  short  account  for  each  term  in 
equation  (35).  The  overriding  factor  in  keeping  AP/P  small  for  the 
terms  on  the  left  hand  side  of  the  equation  is  the  mach  number  M.  As 
long  as  M  <  0.1,  all  of  the  terms  on  the  left  hand  side  give  rise  to  a 
small  value  of  AP/P.  The  viscous  stress  term  on  the  right  hand  side  of 
equation  (35)  can  be  written  as  [29] 


V  •  T  =  pV2U 


'l  3 

r2^ 

•% 

r23r 

3r 

(55) 


where 

U  =  Rr/R. 

Substitution  of  (56)  into  (55)  gives 


(56) 


(57) 
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Prosperetti  [10]  has  shown  that  the  contribution  of  equation  (57)  to  the 

pressure  gradient  in  the  bubble  is  approximately  2yR/PR.  Assuming  that 

-4 

the  viscosity  of  air  at  room  temperature  is  approximately  1.8X10 
Poise,  a  bubble  radius  of  50  microns,  and  a  pressure  of  1  atmosphere, 
one  obtains  a  value  for  AP/P  of  approximately  0.001. 

The  small  value  of  the  viscosity  of  air  also  allows  us  to  neglect 
the  viscous  heating  term  in  the  the  energy  equation  (40).  The  form  of 
this  term  is  given  in  reference  [25]  as 


T 


;pr  =  yR/R. 
rr  3r 


(58) 


The  internal  gas  has  been  modeled  as  an  ideal  gas  in  several  places 
in  the  theory.  This  seems  reasonable  for  the  temperature  and  pressure 
ranges  that  are  of  interest.  This  choice  is  easy  to  justify  since  the 
largest  deviations  from  room  temperature  and  standard  pressure  are 
toward  higher  values  as  the  bubble  collapses.  Callen  [45]  shows  that  at 
higher  temperatures  air  behaves  more  like  that  of  an  ideal  gas.  Because 
of  their  complexity,  more  exact  equations  of  state  would  not  be  easy  to 
incorporate  into  the  theory. 

Finally  we  need  to  discuss  the  error  introduced  when  the  boundary 
condition  (49)  is  used  for  the  energy  equation  (48).  As  long  as  the 
vapor  pressure  of  the  liquid  is  small  compared  to  that  of  the  gas 
pressure,  one  can  assume  that  very  little  vaporization  and  condensation 
takes  place  over  one  cycle.  This  restriction  implies  that  there  is  very 
little  heat  exchanged  on  the  boundary  in  the  form  of  latent  heat.  The 
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other  process  by  which  heat  is  exchanged  on  the  boundary  is  by 

conduction  across  the  interface.  Prosperetti  [10]  has  shown  that  if  T  , 

c 

T  ,  Tm,  and  K  are  the  center  temperature,  surface  temperature,  ambient 
s 

liquid  temperature,  and  the  thermal  conductivity,  then  if  the  heat 
fluxes  in  each  region  across  the  interface  are  set  equal  we  get 

i, 

'2 

’  (59) 

where  1  and  g  denote  the  liquid  and  gas  regions  respectively.  A  typical 
value  for  the  right  hand  side  of  (59)  is  0.001.  This  implies  that  Tg  is 
approximately  equal  to  as  was  assumed  in  the  boundary  condition  (49). 

When  the  liquid  temperature  is  high  enough  so  that  evaporation  and 
condensation  are  important,  one  must  be  careful  in  applying  this 
approximation.  Prosperetti  [10]  gives  a  thorough  discussion  of  this 
problem.  For  the  liquid  temperatures  that  are  of  interest,  however, 
this  does  not  present  a  problem. 

Becc-ase  of  the  complexity  of  the  conservation  equations  and  the 
many  approximations  that  have  been  made,  one  might  ask  if  it  is  really 
worth  all  the  trouble  to  solve  for  P^  in  this  manner.  It  will  become 
evident  in  chapter  3  that  the  extra  work  involved  here  pays  handsome 
dividends  in  terms  of  a  theoretical  model  which  more  closely 
approximates  experimental  observations.  In  chapter  3  the  theories  will 
be  compared  to  each  other  as  well  as  to  experimental  results. 
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Chapter  2 


Numerical  solution  of  the  equations 

A.  Dimensionless  form  of  the  equations 

Perhaps  the  first  step  in  solving  a  system  of  differential 
equations  numerically  is  finding  a  dimensionless  form  of  the  equations 
[46].  In  many  cases  a  dimensionless  form  of  the  equations  is  much  more 
stable  from  a  numerical  point  of  view.  It  is  also  much  easier  to 
compare  the  accuracy  of  different  numerical  methods  if  the  equations  are 
in  a  dimensionless  form.  Before  describing  the  nondimensionalization 
procedure,  let  us  first  recall  the  four  equations  in  the  new  system  with 
the  change  of  variables  given  by  equations  (50)  and  (51) 


fi-il 

••  3*9 

RR  +  -|R2 

r  • 

i  -  r 

_ 

f  d' 

1  +  - 

P(R.t)  -  Ps(t+R/c)  (  R  dP(R.t) 

i  cJ 

3c 

l  CJ 

P  pc  dt  ? 

Y-l  ' 

3t  3t 

Iypr2J 

.3y  "  y3y 

DP 


(52) 


fill 

N 

l  R 

M 

(60) 


and  U  =  dR/dt.  (61) 

In  the  definitions  which  follow,  dimensionless  quantities  will  be 
denoted  with  a  subscript  *.  A  suitable  value  for  a  reference  length  is 
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the  equilibrium  radius  Rq.  The  internal  pressure  at  equilibrium,  P  , 
which  is  related  to  Rq  through  equation  (34),  is  used  as  the  reference 
pressure.  The  reference  time  used  is  the  reciprocal  of  the  driving 
frequency  w  and  the  reference  temperature  used  is  the  ambient  liquid 
temperature  Too*  These  and  other  derived  quantities  are  listed  below  for 
later  reference 


R  ~  V*’  R  ~  V** 

P~  -  V-*»  PS  =  Vs*’  T»  =  T.T*» 

c  =  mRoc*,  t  =  DoPota,  D  *  DoDa. 


(62) 


The  appropriate  form  of  the  different  derivatives  in  our  system  of 
equations  is  calculated  next.  The  time  derivatives  of  the  radius  are 


R  *  and  R  =  u)2RqRa  ,  (63) 

where  the  dots  over  starred  variables  denote  differentiation  with 
respect  to  the  dimensionless  time  t*.  The  other  derivatives  necessary 
to  complete  the  non-dimensionalization  process  are 
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where  several  steps  have  been  omitted,  particularly  in  the  last 
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of  the  form 
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K(T)  =  5.528T  +  1165.  (68) 

Equation  (68)  was  determined  by  a  simple  linear  regression  routine  on 
the  experimental  data  given  by  Weast  [47]  (see  figure  1).  This  linear 
equation  for  K(T)  enables  one  to  integrate  equation  (50)  analytically 
and  obtain  a  value  for  T  from  the  resulting  quadratic  equation. 

This  concludes  the  nondimensionalization  process  for  the  system  of 
equations.  The  details  on  how  to  solve  these  equations  is  the  subject 
of  the  next  section  of  this  chapter. 


20000 


degrees  kelvin. 
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B.  Finite  difference  equations 

The  system  of  dimensionless  equations  derived  in  section  A  of  this 
chapter  is  now  solved  using  numerical  methods.  Because  of  the 
complexity  of  the  system  and  the  coupling  between  the  equations,  a 
multistep  method  is  used.  These  types  of  methods  offer  good  accuracy 
without  being  as  time  consuming  as  some  of  the  single  step  algorithms 
such  as  the  Runge-Kutta  method. 

The  algorithm  used  is  as  follows.  A  simple  predictor-corrector 
multistep  method  is  applied  to  our  system  of  equations  [50].  This 
technique  can  be  illustrated  by  letting 

dX 

—  =F(X).  (69) 

The  vector  _X  denotes  the  quantities  radius,  velocity,  pressure,  and 
temperature.  FOO  is  the  functional  form  of  the  time  derivatives  of  the 
aforementioned  quantities.  Using  finite  differences  one  can  write 
equation  (69)  as 

Xn+1  =  Xn  +  -y-[F(Xn)  +  F(Xn+1)].  (70) 

However,  one  cannot  solve  this  equation  explicitly  for  _Xn+^  because  of 
the  complexity  of  the  system  of  equations.  The  predictor-corrector 
scheme  overcomes  this  problem  however. 

In  the  predictor  step  an  approximate  value  for  _X°+^  is  computed 
using  an  Euler  method  [48]  of  the  form 

Xn+1  =  XU  +  AtF(Xn). 


(72) 
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This  predicted  value  of  _X  is  used  in  equation  (70)  to  get  a  more 
accurate  value  of  _Xn+*.  In  actual  practice  this  corrector  step  is 
repeated  a  second  time  for  additional  accuracy.  One  could  repeat  this 
process  over  and  over,  but  experience  has  shown  that  two  iterations  is 
sufficient  for  our  needs. 

If  one  finds  certain  initial  conditions  and  parameter  settings 
which  cause  the  corrector  process  to  converge  slowly,  a  more  advanced 
predictor-corrector  algorithm  such  as  the  Adams-Moulton  method  [48]  can 
be  used.  For  the  conditions  that  are  of  interest,  however,  the  simple 
method  suffices. 

Although  the  techniques  outlined  above  are  straightforward  in 
principle,  they  are  rather  tedious  to  apply  because  of  the  complexity  of 
the  system  of  equations.  In  the  remainder  of  this  section,  a  step  by 
step  procedure  on  how  to  implement  this  method  for  this  system  of 
equations  is  discussed.  Those  not  interested  in  this  detailed  guide  to 
the  solution  of  the  equations  are  encouraged  to  skip  to  the  flow  chart 
at  the  end  of  this  section  without  loss  of  continuity. 

The  first  phase  of  the  method  is  to  predict  new  values  of  the 
variables  using  the  Euler  equation  (71).  For  convenience  the  subscript 
*  notation  for  dimensionless  equations  (66)  is  dropped.  Predicted 
values  will  be  denoted  with  a  tilda  over  them.  The  first  step  is  finding 


a  predicted  pressure  at  the  n+1  time  step. 
dP/dt  at  the  time  step  by 

=  _3_  (y- 1 ) »,  f0~V 

Rn  Rn  l  . 


This  requires  the  evaluation 
-  yPnUn) ,  (72) 


Chapter  3 


Computational  results 

A.  Radius  verses  time  curves 

Of  all  the  information  obtained  from  the  computer  programs,  perhaps 
the  simplest  to  understand  and  analyze  is  the  radius  verses  time  data. 
This  information  gives  a  numerical  picture  of  the  bubble  at  each  time 
step.  One  can  compare  these  radius-time  curves  for  different  models  at 
a  variety  of  initial  conditions  and  make  some  conclusions  as  to  where 
the  models  agree  and  disagree.  Because  some  models  such  as  the  exact 
formulation  require  large  amounts  of  computer  time,  it  is  useful  to  know 
when  a  simpler  model  gives  results  that  are  sufficiently  accurate  for 
one's  needs.  For  this  reason  we  not  only  show  results  from  the  new 
formulation  but  also  include  the  results  from  other  methods  as  well. 
One  can  also  obtain  information  about  the  phase  difference  between  the 
driving  pressure  and  the  bubble  oscillations.  This  phase  difference  is 
related  to  the  damping  of  the  system  in  a  complicated  manner.  Finally 
one  can  find  the  spectrum  of  the  radius-time  curves  by  a  simple  fourier 
analysis  of  the  data.  This  fourier  analysis  of  the  curve  shows  which 
harmonic  modes  are  dominant  in  the  bubble's  radiated  noise  signal. 

In  figure  5,  the  radius  verses  time  curves  are  compared  for  an 
equilibrium  radius  of  50  microns,  a  driving  pressure  amplitude  of  0.4 
bar,  and  a  driving  frequency  of  0.4  times  the  linear  resonance 
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temperature  of  bubbles  at  known  driving  pressures  and  equilibrium  radii. 
An  investigation  of  this  nature  is  presently  under  way  and  the  results 
should  be  available  shortly. 

In  closing  we  should  mention  that  experience  has  shown  that  by 
examining  a  few  cycles  of  the  numerical  solution,  one  can  easily  decide 
whether  equation  (70)  is  stable  or  unstable.  For  the  stable  solutions 
at  high  amplitudes,  one  must  wait  for  new  experimental  results  before 
deciding  how  well  the  new  theory  models  a  real  cavitation  bubble  in  this 
extreme  region. 
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D.  Applicability 

In  the  next  chapter  the  values  given  by  the  new  equations  are  shown 
to  be  more  accurate  than  the  ones  given  by  the  polytropic  approximation. 
This  is  shown  by  examining  the  predicted  levitation  numbers  of  each 
method  and  comparing  them  to  experimental  ones  obtained  by  Crum  and 
Prosperetti  [  19]  . 

Unfortunately,  the  levitation  number  data  is  still  at  rather  low 
amplitude  and  one  cannot  tell  how  accurate  the  solutions  are  in  high 
amplitude  regions.  One  possible  experimental  method  to  test  the 
validity  of  the  numerical  solution  at  high  amplitudes  is  outlined  below. 
One  could  use  the  light-scattering  technique  of  Hansen  [51]  with  a  very 
fast  photo-diode  and  a  good  low  noise  amplifier  attached.  This  would 
produce  a  signal  that  is  related  to  the  radius  of  the  bubble  by  the 
relations  given  by  Hansen.  Once  properly  converted  to  a  radius  verses 
time  signal  one  could  compare  the  numerical  solution  at  high  amplitudes 
to  this  experimental  one.  Work  on  this  experiment  is  presently  under 
way  and  when  completed  it  should  provide  a  definitive  test  for  all 
nonlinear  bubble  oscillation  equations. 

Another  test  for  large  amplitude  oscillations  is  determining  the 
interior  bubble  temperatures  experimentally  by  observing  sono- 
luminescence  and  comparing  them  to  the  numerical  values.  Dissociation 
of  water  molecules  into  hydroxyl  free  radicals  and  hydrogen  starts 
taking  place  at  well  defined  temperatures  [7],  which  should  enable 
experimentalists  to  obtain  very  accurate  data  on  the  interior 
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figure  J.  Normalized  bubble  radius  as  a  function  or  time 
for  the  analytic  first  order  solution  and  the  new  exact 
formulation.  R  =50  microns,  Pa=0.01  bar,  and  (0=10  /2. 
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insight  as  to  how  well  the  numerical  solution  approximates  an  exact 
analytic  solution  to  the  system  of  equations  if  one  were  available.  One 
can  make  several  tests,  however,  to  increase  one's  confidence  in  the 
method  used.  As  a  preliminary  test  the  numerical  solution  for  the 
radius  as  a  function  of  time  is  compared  to  an  analytical  solution  which 
uses  the  polytropic  approximation  and  is  valid  to  first  order  [15].  Of 
course  this  comparison  can  only  be  made  at  low  driving  amplitudes  since 
the  analytical  solution  is  only  accurate  in  the  low  amplitude  limit. 
Figure  3  shows  this  comparison  for  a  50  micron  bubble  driven  at  half  its 
linear  resonance  frequency  at  a  pressure  amplitude  of  0.01  bar.  The 
figure  indicates  that  the  numerical  solution  is  nearly  identical  to  the 
analytical  one  when  driving  amplitudes  are  small.  This  reassures  us 
that  a  fundamental  mistake  has  not  been  made  in  either  deriving  the 
system  of  equations  or  in  solving  them  numerically. 

For  more  moderate  amplitudes,  a  comparison  is  made  between  the  new 
results  and  numerical  solutions  of  the  Rayleigh-Plesset  equation  with  a 
polytropic  approximation,  using  a  fourth  order  Runge-Kutta  method.  The 
program  used  to  do  this  calculation  was  first  written  by  A.  Prosperetti. 
An  updated  version  with  modifications  made  by  the  author  is  included  in 
appendix  B.  Figure  4  shows  a  comparison  between  the  two  models  for  a  50 
micron  bubble  driven  at  half  its  linear  resonance  frequency  and  a 
pressure  of  0.50  bar.  For  this  case  there  is  a  small  difference  between 
the  two  solutions,  however,  they  are  similiar  enough  to  reassure  us  that 
the  new  solution  is  working  properly  at  this  amplitude. 


C.  Error  analysis 

A  formal  error  analysis  of  the  complete  system  of  equations  in  all 

its  detail  would  be  quite  difficult.  One  can,  however,  give  a  more 

general  error  analysis  for  the  basic  methods  used,  that  is,  for 

equations  (70)  and  (71)  and  the  Crank-Nicolson  implicit  method. 

The  Euler  method  or  equation  (71),  which  is  used  as  the  predictor 

2 

in  the  algorithm,  has  a  local  or  truncation  error  of  order  (At)  and  a 
global  or  cumulative  error  of  order  (At).  Equation  (70)  is  known  as  the 
Modified  Euler  method  or  the  Euler  predictor-corrector  method.  The 

3 

truncation  error  of  the  Modified  Euler  method  is  of  order  (At)  and  the 

2 

cumulative  error  is  of  order  (At)  .  For  completeness  one  should  also 
mention  the  Adams-Moulton  method  referred  to  in  the  previous  section. 
This  method  has  a  truncation  error  of  order  (At)'’  and  a  cumulative  error 

4 

of  order  (At)  ,  which  is  the  same  as  a  fourth  order  Runge-Kutta  method. 

The  Crank-Nicolson  method  is  embedded  in  equation  (70).  This 

method  is  an  implicit  one  and  because  of  the  judicious  choice  of  the 

constants  multiplying  the  finite  difference  approximations,  the  method 

2  2 

has  a  truncation  error  of  order  (At)  and  (Ay)  •  Since  this  method  is 
always  stable,  one  knows  that  errors  created  at  each  time  step  must 
decay  exponentially  [48].  One  of  the  nice  features  of  damped  driven 
nonlinear  oscillators  is  that  as  long  as  the  algorithms  are  stable  the 
numerical  solution  is  usually  being  forced  toward  the  exact  solution  of 
the  equations. 

The  error  analysis  above  does  not,  however,  give  any  physical 
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(88) 

Finally  the  corrected  values  for  the  dimensional  temperature  and  the 
thermal  diffusivity  are  found  by  replacing  the  tilda  variables  in 
equations  (77)  and  (79)  with  the  corrected  variables  denoted  by  carets. 

In  order  to  increase  the  accuracy  at  the  new  time  step  even  more, 
a  second  correction  exactly  the  same  as  the  one  just  described  is  made. 
It  is  these  second  corrections  that  are  used  for  the  values  at  the  n+1 
time  step.  At  this  point  in  the  computer  program,  all  of  the  variables 
that  are  of  interest  are  written  out  and  then  the  variable  names  are 
updated.  The  entire  process  is  then  repeated  starting  at  equation  (72). 

A  listing  of  the  computer  program  is  provided  in  appendix  A. 
Figure  2  is  a  flowchart  for  the  program  and  provides  a  qualitative  view 
of  how  the  system  of  equations  is  solved.  The  error  analysis  for  the 
numerical  methods  used  is  found  in  the  next  section  of  this  chapter. 
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and  the  corrected  pressure  derivative  from  equation  (66c)  is  given  by 
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Finally  the  corrected  values  for  the  dimensional  temperature  and  the 
thermal  diffusivity  are  found  by  replacing  the  tilda  variables  in 
equations  (77)  and  (79)  with  the  corrected  variables  denoted  by  carets. 
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that  are  of  interest  are  written  out  and  then  the  variable  names  are 
updated.  The  entire  process  is  then  repeated  starting  at  equation  (72). 

A  listing  of  the  computer  program  is  provided  in  appendix  A. 
Figure  2  is  a  flowchart  for  the  program  and  provides  a  qualitative  view 
of  how  the  system  of  equations  is  solved.  The  error  analysis  for  the 
numerical  methods  used  is  found  in  the  next  section  of  this  chapter. 
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The  next  step  in  the  corrector  method  is  finding  the  corrected 
bubble  wall  velocity.  Again,  equation  (70)  is  applied  as  before.  One 
must  be  careful,  however,  because  the  F^(JCn+^)  term  has  both  predicted 
and  corrected  velocities  in  it.  The  appropriate  form  of  equation  (70), 


before  solving  for  the  corrected  velocity  explicitly,  is  given  by 
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After  equation  (85)  is  solved  for  the  corrected  velocity  explicitly,  the 
corrected  radius  at  the  new  time  step  is  found  by  applying  equation  (70) 


which  gives 


-n+1  ,  Rn  +  M(un+0"+1). 


Using  the  corrected  values  of  temperature,  velocity,  and  radius 
one  finds  that  the  corrected  pressure  from  equation  (70)  is  given  by 
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variable,  one  starts  the  corrector  algorithm  using  the  method  of 
equation  (70).  The  corrected  variables  will  be  denoted  by  a  caret  above 
them.  The  corrector  method  begins  by  solving  for  a  corrected 
temperature  array  at  the  n+1  time  step.  This  time,  however,  it  is  not  a 
simple  matter  to  solve  for  the  new  values  because  one  must  use  a  Crank- 
Nicolson  method  {48]  to  advance  the  temperature  array  in  time.  This  is 
necessary  because  equation  (70)  is  accurate  to  second  order  and  each  of 
the  embedded  algorithms  must  be  at  least  this  accurate  or  one  has 
defeated  the  purpose  of  using  equation  (70).  The  Crank-Nicolson  method 
is  an  implicit  method  which  requires  the  solution  of  a  tridiagonal 
system  of  linear  equations  for  the  corrected  temperature  array.  In 
order  to  see  this  more  clearly,  the  energy  equation  is  written  in  a  form 
with  the  unknowns  and  their  coefficients  on  the  left  hand  side  and  all 
known  quantities  on  the  right  hand  side.  As  before,  the  first  node  has 
a  special  form  given  by 
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The  general  expression  for  the  remaining  nodes  from  i  =*  2  to  N  is  given 
by  equation  (84)  below.  To  solve-  his  tridiagonal  system  for  the 
corrected  temperatures,  an  algorithm  by  Forsythe,  Malcolm,  and  Moler 
[49]  was  modified  for  speed.  This  algorithm  is  a  subroutine  called 


TRIDG  in  our  main  computer  program  which  is  listed  in  appendix  A. 
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thermal  diffusivity  D  at  each  node.  This  requires  that  the  dimensional 
temperature  be  found  first  by  the  method  mentioned  in  section  A  of  this 
chapter.  For  each  node  the  dimensional  temperature  is  related  to  the 
dimensionless  temperature  by 
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where 

a  =  5.528^/^). 
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The  predicted  diffusivity  at  node  i  is  given  by 

5”+1  +  (ai"+1  +  1165/K(TJ)T"+i/Pn+1.  (79) 

We  can  now  return  to  finding  new  values  for  the  main  variables.  A 
predicted  value  for  the  radius  is  obtained  from 
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To  find  a  predicted  velocity,  one  first  solves  equation  (66a)  for  dU/dt 
and  evaluates  it  at  the  nC^  time  step.  After  evaluating  this 
derivative,  the  predicted  velocity  is  found  from  the  Euler  formula  by 
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The  predictor  section  is  concluded  once  a  new  value  of  dP/dt  is  found  by 
•n+1  _  3 
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Having  found  a  predicted  value  at  the  new  time  step  for  each 
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where  subscripts  denote  the  node  or  spatial  position  and  superscripts 
denote  the  time  step.  Using  Euler’s  method  one  finds  the  predicted 
pressure  is  given  by 


pn+1 
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The  next  step  involves  the  temperature  array  in  the  bubble.  First 
one  has  to  find  the  temperature  derivatives  with  respect  to  time  at  each 
node  in  the  bubble.  The  r=0  node  is  a  special  case  since  the  Laplacian 
operator  is  different  in  the  limit  as  r  approaches  0  than  at  any  nonzero 
value  of  r  [46].  The  center  node  or  i=l  node  is  found  by 

in 
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The  general  equation  for  the  remaining  nodes  from  i=2  to  N  (  N+l  is  the 
bubble  wall  which  was  assumed  to  be  at  constant  ambient  temperature  )  is 
given  by 
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Again  using  the  Euler  method  one  finds  predicted  temperatures  at  the  new 
time  step  n+l  by 
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At  this  point  it  is  appropriate  to  find  predicted  values  of  the 
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frequency.  The  Rayleigh-Plesset  Polytropic  equation  is  denoted  by  RPP, 
the  Rayleigh-Plesset  Exact  pressure  formulation  by  RPE,  and  the 
Compressible  Ebcact  pressure  system  of  equations  by  CE.  The  heavy  curve 
in  figure  5  represents  the  RPE  and  CE  methods  which  are  practicably 
identical.  The  compressible  equation  is  only  necessary  when  the  driving 
pressure  amplitude  is  very  large  and  collapse  velocities  are  very  high. 
It  takes  only  a  small  amount  of  extra  computer  time  to  use  the 
compressible  equation  in  the  exact  formulation,  so  we  will  use  the  CE 
method  and  omit  from  further  discussion  the  RPE  method. 

Next,  a  comparison  is  made  between  the  solutions  given  by  the  CE 
and  RPP  methods  for  a  variety  of  bubble  sizes,  driving  pressure 
amplitudes,  and  driving  frequencies.  Figure  6  shows  results  for  an 
equilibrium  radius  of  50  microns,  a  driving  pressure  amplitude  of  0.4 
bar  and  a  driving  frequency  of  0.8  times  the  linear  resonance  frequency. 
The  CE  method  predicts  a  larger  maximum  value  of  R/Rq  than  does  the  RPP 
method.  In  figure  7  the  equilibrium  radius  is  50  microns,  the  driving 
pressure  amplitude  is  0.6  bar,  and  the  driving  frequency  is  0.8  times 
the  linear  resonance  frequency.  Again  the  CE  method  predicts  a  larger 
maximum  value  of  R/Rq  than  does  the  RPP  method.  Note  also  the  phase 
difference  between  the  two  curves  which  is  easily  measurable.  These 
curves  have  very  large  amplitudes,  however,  and  the  limits  of  validity 
of  the  equations  may  have  been  reached  or  exceeded.  One  can  only  wait 
for  new  experimental  data  to  determine  just  how  large  an  amplitude  can 
be  used  before  appreciable  error  occurs  in  the  results. 
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Figure  5.  Normalized  bubble  radius  as  a  function  of  time 
for  the  RPP,  RPE,  and  CE  methods.  RPE  and  CE  above  the  RPP 
curve.  R  =50  microns,  Pa=0.4  bar,  and  di=0.4t,i  . 
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Normalized  bubble  radius  as  a  function  of  time 
'  and  CE  methods.  CE  above  RPP.  R  =50  microns 
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Next,  several  figures  for  10  micron  bubbles  are  examined.  Figure  8 
represents  a  10  micron  bubble  with  a  driving  pressure  amplitude  of  0.6 
bar  and  a  driving  frequency  of  0.4  times  the  linear  resonance  frequency. 
Again  the  pattern  is  the  same  as  before,  with  the  CE  method  showing  less 
damping  than  the  RPP  method.  One  can  also  detect  a  difference  in  the 
magnitude  of  the  second  harmonic  resonance  by  inspection  of  the  figure 
(  the  small  peaks  between  the  larger  peaks  are  mainly  due  to  the  second 
harmonic  ).  Figure  9  concludes  this  set  of  results  with  an  equilibrium 
radius  of  10  microns,  a  driving  pressure  of  0.6  bar,  and  a  driving 
frequency  of  0.8  times  the  linear  resonance  frequency.  Again  one  gets  a 
larger  response  in  the  amplitude  of  the  pulsation  because  the  driving 
frequency  is  in  the  vicinity  of  the  main  resonance,  whereas  the  curves 
at  the  0.4  frequency  are  near  the  second  harmonic.  In  the  next  section 
it  is  shown  that  the  response  due  to  these  harmonic  resonances  is  much 
less  than  that  at  the  main  resonance. 

These  figures  represent  only  a  few  of  the  unlimited  number  of 
possible  combinations  of  initial  conditions.  They  do,  however,  give  an 
idea  of  the  differences  in  the  two  theories  presented.  With  careful 
observation  one  can  detect  differences  in  the  pulsation  amplitude,  the 
phase  shifts  between  the  driving  pressure  and  the  bubble  response,  and 
in  the  magnitude  of  each  harmonic  if  a  fourier  analysis  is  carried  out 
for  the  curves  generated. 
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Figure  9.  Normalized  bubble  radius  as  a  function  of  time 
for  the  RPP  and  CE  methods.  CE  above  RPP.  R  =10  microns 
Pa=0.6  bar,  and  u)=0.8u)  .  ° 
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B.  Frequency  response  curves 

One  of  the  best  ways  of  showing  how  a  bubble  responds  to  a  driving 
force  under  a  variety  of  initial  conditions  and  parameter  settings  is 
through  the  use  of  frequency  response  curves  such  as  those  given  by 
Lauterborn  [18,52]  anj  Cramer  [17].  These  graphs  show  the  nonlinear 
resonance  peaks  which  help  to  characterize  a  bubble.  This  information 
is  of  interest  to  those  examining  the  biological  effects  of  ultrasound 
since  the  onset  of  subharmonic  frequencies  is  often  used  as  an 
indication  of  the  presence  of  violent  cavitation.  Therefore,  it  is 
important  to  know  precisely  at  what  fraction  of  the  bubble’s  resonance 
frequency  that  these  different  modes  of  oscillation  prevail  for 
different  driving  pressure  amplitudes. 

In  the  figures  that  follow,  results  are  shown  for  a  variety  of 
conditions  using  the  new  theory.  These  results  are  compared  to  those  by 
Lauterborn  and  Cramer  to  determine  where  the  two  theories  agree  and 
disagree.  Figures  10-12  are  frequency  response  curves  generated  by  the 
new  theory  for  normalized  frequencies  between  0.1  and  1.0.  The  step 
size  used  in  the  frequency  domain  is  0.01  times  the  resonance  frequency 
of  that  particular  bubble.  Figure  10  shows  results  for  a  10  micron 
bubble  driven  at  0.3,  0.5,  and  0.7  bar.  Figure  11  shows  results  for  a 
50  micron  bubble  driven  at  0.2,  0.3,  0.4,  and  0.5  bar.  Finally  figure 
12  shows  results  for  a  100  micron  bubble  driven  from  0.1  to  0.7  bar  in 
0.1  bar  intervals. 

Each  of  the  fourteen  curves  described  above  took  approximately  one 
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day  of  cpu  time  on  a  Digital  PDP  11/73  microcomputer.  Of  course  the 
curves  which  show  a  large  amplitude  response  took  longer  than  those  with 
small  responses  because  of  the  difference  in  step  sizes  used  in  the  time 
and  space  domains.  This  requirement  of  large  amounts  of  computer  time 
has  somewhat  limited  the  scope  of  this  investigation;  this  point  will  be 
addressed  again  in  chapter  4. 

In  figures  13  and  14  new  results  are  compared  with  those  by 

Lauterborn  and  Cramer.  Figure  13  shows  the  frequency  response  curve  for 

a  10  micron  bubble  driven  at  0.7  bar.  The  asterisks  represent  points 

taken  from  the  peaks  of  a  figure  by  Lauterborn  [52]  (  the  value  at 

f/f  *1.0  is  included  as  well  ).  Note  that  there  is  a  shift  in  both 
o 

amplitude  and  frequency  between  the  two  theories.  As  f/fQ  decreases  the 
two  theories  give  results  that  are  quite  close  to  one  another.  Figure 
14  shows  the  new  frequency  response  curve  for  a  100  micron  bubble  driven 
at  0.7  bar.  Here  the  asterisks  represent  points  taken  from  the  peaks  of 
a  figure  by  Cramer  [17].  Again,  one  sees  the  same  trend  as  before,  with 
small  values  of  f /f q  giving  even  closer  results  than  with  the  10  micron 
bubble.  These  figures  indicate  that  the  new  theory  predicts  less 
damping  in  the  region  around  the  main  resonance  peak  and  some  of  the 
subharmonic  peaks  as  well.  From  the  comparison  to  levitation  number 
data  in  the  next  section  one  is  inclined  to  believe  that  this  new  theory 
is  more  accurate.  The  definitive  test,  however,  is  experimental  data 
which  is  not  yet  available. 

Close  examination  of  figure  13  and  Lauterborn's  figure  3  [52] 
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normalized  frequency  for  the  CE  method.  R  =10  microns, 
0.7  bar,  asterisks  denote  resonance  peaks  from  a  figure 
Lauterborn  [ 52 } . 
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the  gas  inside  starts  to  increase  due  to  heat  conduction  across  the 
bubble  wall.  For  small  bubbles,  the  value  of  the  radius  where  the 
minimum  temperature  is  reached  is  near  Rq,  and  for  large  bubbles  it  is 
near  the  maximum  radius  attained. 

The  figures  above  indicate  that  under  the  proper  conditions, 
cavitating  bubbles  can  produce  temperatures  high  enough  for  free  radical 
formation.  This  should  cause  some  concern  for  those  using  ultrasonic 
instruments  for  medical  purposes.  One  way  to  avoid  this  problem  is  to 
stay  away  from  driving  pressures  and  frequencies  that  cause  large 
pulsation  amplitudes.  To  avoid  these  regions,  one  must  first  know  the 
size  of  the  bubbles  in  the  system  of  interest.  This,  however,  is  a 
difficult  problem  in  itself  and  will  not  be  addressed  further. 

Next,  the  center  temperature  of  a  bubble  is  examined  as  a  function 
of  time  to  see  how  long  it  maintains  these  high  temperatures.  Figures 
22-24  show  center  temperatures,  as  a  function  of  time  for  1,  10,  and  100 
micron  bubbles  driven  at  0.4  and  0.8  times  the  linear  resonance 
frequency  of  the  particular  bubble.  Each  figure  indicates  that  high 
temperatures  occur  only  for  a  short  period  of  time  while  the  bubble  is 
near  its  minimum  size.  Figure  22  shows  a  1  micron  bubble  driven  at  1.4 
bar,  figure  23  a  10  micron  bubble  driven  at  0.6  bar,  and  figure  24  a  100 
micron  bubble  driven  at  0.45  bar.  The  pulsation  amplitude  for  the  10 
micron  bubble  is  much  smaller  than  either  the  l  or  100  micron  bubble. 
Thus,  the  maximum  temperature  attained  by  the  10  micron  bubble  is  much 
less  than  in  either  of  the  other  two  bubbles. 
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Values  for  the  minimum  and  maximum  radius  are  shown. 
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volume  and  time,  using  the  temperature  profiles  and  an  appropriate 
ionization  probability  function  as  the  integrand. 

Figures  19-21  are  temperature  profiles  for  1,  10,  and  100  micron 
bubbles  at  different  stages  of  oscillation.  The  minimum  and  maximum 
value  of  the  radius  for  each  bubble  during  one  cycle  is  chosen  as  the 
point  where  the  temperature  profiles  are  taken.  The  driving  frequency 
for  each  figure  is  0.4  times  the  linear  resonance  frequency  of  the 
bubble.  The  temperatures  shown  are  normalized  with  respect  to  the 
ambient  temperature  of  293  degrees  kelvin.  Figure  19  shows  a  1  micron 
bubble  driven  at  2.1  bar.  Note  that  when  the  radius  is  a  minimum,  the 
center  temperature  is  approximately  3.2  times  the  equilibrium 
temperature.  Figure  20  shows  a  10  micron  bubble  driven  at  0.9  bar. 
Here  the  maximum  temperature  attained  is  approximately  4.9  times  the 
equilibrium  temperature.  Finally,  figure  21  shows  a  100  micron  bubble 
driven  at  0.8  bar.  For  this  case  the  maximum  temperature  attained  is 
approximately  8.1  times  the  equilibrium  value. 

In  each  of  these  figures  the  maximum  value  of  R/Rq  Is  approximately 
2.0.  Thus,  the  small  bubbles  behave  much  more  isotherraally  than  do  the 
large  ones.  Since  the  bubble  wall  moves  slowly  at  large  values  of  R/Rq, 
the  interior  temperature  does  not  drop  to  the  level  that  might  be 
expected  at  the  maximum  value  of  the  radius.  Instead,  a  minimum 
temperature  is  reached  when  the  radius  of  the  bubble  is  somewhere 
between  its  equilibrium  value  and  maximum  value  during  the  expansion 
cycle.  As  the  bubble  approaches  its  maximum  size,  the  temperature  of 
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D.  Thermodynamics  of  the  bubble  interior 

As  has  been  mentioned  previously,  the  internal  pressure  and 
temperature  of  a  bubble  are  of  particular  interest  to  those  studying  the 
physiological  effects  of  acoustic  cavitation.  There  are  several 
mechanisms  by  which  a  cavitating  bubble  can  damage  its  surroundings. 
The  best  known  of  these  mechanisms  is  that  of  destruction  by  the  bubble 
wall  itself.  This  area  has  been  thoroughly  studied  because  of 
cavitation  damage  to  propellers  on  ships.  However,  we  will  not  address 
this  method  of  cavitation  damage  here.  More  recently,  experimentalists 
have  found  that  temperatures  inside  the  bubble  can  become  high  enough  to 
produce  free  radicals.  The  free  radicals  produced  at  these  high 
temperatures  can  pose  a  serious  threat  to  biological  systems.  One 
method  of  observing  this  phenomena  is  to  look  for  sonoluminescence  from 
cavitating  bubbles. 

The  exact  formulation  for  the  internal  pressure  of  the  bubble  has 
within  itself  a  temperature  array  for  the  bubble  interior.  Thus,  one 
solves  for  the  internal  temperature  of  a  bubble  as  a  function  of  both 
position  and  time.  It  is  important  to  know  how  many  free  radicals  are 
produced  in  a  bubble  over  one  cycle  for  a  given  set  of  parameters.  One 
can  estimate  this  number  from  the  temperature  profiles  the  computer 
program  generates.  One  needs  to  know  what  fraction  of  the  bubble's 
volume  reaches  the  critical  temperature  for  free  radical  formation  and 
how  long  it  remains  at  or  above  this  temperature.  A  theoretical  number 
of  free  radicals  could  be  obtained  by  numerically  integrating  over 
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Figure  18.  Levitation  number  as  a  function  of  the  normalized 
radius  for  the  CE  and  RPP  methods.  Freq . =22 . 2kHz ,  Pa=0.190 
bar,  asterisks  denote  experimental  data  from  Crum  [19]. 
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Figure  17.  Levitation  number  as  a  function  of  the  normalized 
radius  for  the  CE  and  RPP  methods.  Freq.=22. 2kHz,  Pa=0.155 
bar,  asterisks  denote  experimental  data  from  Crum  [19]. 
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Figure  16.  Levitation  number  as  a  function  of  the  normalized 
radius  for  the  CE  and  RPP  methods.  Freq.=22.2kHz,  Pa=0.095 
bar,  asterisks  denote  experimental  data  from  Crum  [19]. 
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as  1%.  When  this  experiment  is  complete  one  should  be  able  to  determine 
the  accuracy  of  the  new  theoretical  model  precisely. 

Figures  16-18  show  results  for  a  frequency  of  22.2  kHz  and  pressure 
amplitudes  of  0.095,  0.155,  and  0.190  bar  respectively.  In  each  figure 
the  asterisks  denote  experimental  points  from  reference  [19].  The  curve 
giving  the  best  fit  to  the  experimental  data  in  each  figure  is  the  new 
theoretical  model.  From  these  figures  one  can  conclude  that  the  new 
theory  is  more  accurate  than  the  old  polytropic  exponent  theory  for  the 
parameter  range  investigated. 


the  acoustic  force  is  balanced  by  an  average  buoyancy  force  given  by 


fb  =  '|TTRoP8<[R(t)/R0l3>' 


(92) 


where  g  is  the  acceleration  due  to  gravity.  Equating  these  two  forces 
gives 


IVP  |<[R(t)/R  ]3Costot>  =  pg<fR(t)/R  ]3>. 
**  o  o 


(93) 


Since  pg  is  the  hydrostatic  pressure  gradient,  one  can  rearrange 

equation  (93)  so  that  the  right  hand  side  is  the  levitation  number  as 

defined  above.  Thus  equation  (93)  becomes 

3 

<[R(t)/R  ]  cosu)t>  pg 

2  -  =  - - -  =  L  .  (94 ) 

<[R(t)/RQ]J>  | VPA| 

The  left  hand  side  of  equation  (94)  can  easily  be  evaluated  for  a 
numerical  solution  to  equations  modeling  bubble  oscillations.  From 
this  equation  one  can  compute  levitation  numbers  for  both  the  CE  method 
and  the  RPP  method  and  compare  the  results  to  the  experimental 
levitation  numbers  obtained  by  Crum.  The  experimental  points  are 
obtained  in  a  rather  straightforward  manner  and  the  technique  is 
discussed  in  detail  in  reference  [19].  There  may  be  as  much  as  a  10% 
error  in  obtaining  the  levitation  numbers  experimentally.  However,  it  is 
assumed  that  the  experimental  data  is  accurate  enough  to  determine  which 
theoretical  model  best  predicts  the  time  evolution  of  an  oscillating 
bubble  when  pulsation  amplitudes  are  not  too  large.  A  more  automated 
experimental  procedure  is  being  devised  which  could  give  results  with 
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C.  Levitation  numbers 

The  ratio  of  the  hydrostatic  pressure  gradient  on  a  bubble  to  the 

acoustic  pressure  gradient  has  been  referred  to  as  the  levitation  number 

L  by  Crum  and  Prosperetti  [19].  Before  using  this  definition  with  the 
e 

new  theory,  let  us  first  recall  where  the  expression  arises  in  the  first 
place. 

The  primary  Bjerknes  force  or  acoustic  radiation  pressure  force  on 
a  bubble  is  given  by  Crum  [53]  as 

V£’c)  “  -<V(t)VP(r,t)>,  (89) 

where  V(t)  is  the  time  dependent  volume  of  the  bubble,  P(jr,t)  is  the 
space  and  time  dependent  pressure  outside  the  bubble,  and  the  angle 
brackets  denote  the  time  average  of  the  quantity  inside.  If  one  uses  a 

levitation  cell  like  the  one  used  by  Crum,  then  the  pressure  PCjr, t)  in 

the  cell  is  of  the  form 

P(£,t)  =  Pm  -  PA(z)coswt,  (90) 

where  P,(z)  is  the  spatially  dependent  amplitude  of  the  acoustic 
A 

standing  wave  in  the  cell.  Substitution  of  equation  (90)  into  (89) 
gives  the  magnitude  of  the  acoustic  force  as 

FA  =  3irRolvpAl<[R(t>/R0]3coswt>,  (91) 

where  Rq  is  the  equilibrium  radius  of  the  bubble. 

When  the  bubble  remains  in  a  fixed  position  in  the  levitation  cell, 
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reveals  missing  ultraharmonic  resonance  peaks  in  figure  13.  For 
example,  Lauterborn's  figure  shows  a  3/2  ultraharmonic,  whereas  figure 
13  does  not.  In  order  to  be  sure  that  these  peaks  are  not  being  skipped 
by  too  large  a  step  size  in  the  frequency  domain,  another  calculation 
is  made  with  twice  as  many  points  in  the  frequency  domain  than  before. 
This  result  is  presented  in  figure  15  for  a  10  micron  bubble  driven  at 
0.7  bar  along  with  some  partial  curves  at  0.8  and  0.9  bar.  Points  taken 
from  ultraharmonic  peaks  in  Lauterborn's  figure  are  denoted  by 
asterisks.  Again  the  ultraharmonic  resonances  are  missing  from  our 
curve  at  0.7  bar,  but  they  begin  to  appear  as  the  pressure  is  increased. 
Qualitatively,  the  two  theories  give  similiar  results,  but  there  are 
significant  quantitative  differences  which  can  be  important  in  certain 
applications  as  mentioned  previously. 

Although  there  are  good  indications  that  the  new  theory  is  more 
accurate  in  some  parameter  domains,  one  should  wait  for  future 
experimental  results  before  deciding  which  theory  best  predicts  the 
frequency  response  of  the  bubble  in  these  high  amplitude  regions. 
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Figure  22.  Center  temperature  as  a  function  of  time  for 
a  1  micron  bubble  driven  at  1.4  bar,  and  frequencies  of 
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Figure  23.  Center  temperature  as  a  function  of  time  for 
a  10  micron  bubble  driven  at  0.6  bar,  and  frequencies  of 
0.4  and  0.8  times  the  linear  resonance  frequency. 
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Figure  24.  Center  temperature  as  a  function  of  time  for 
a  100  micron  bubble  driven  at  0.45  bar,  and  frequencies 
of  0.4  and  0.8  times  the  linear  resonance  frequency. 


To  complete  this  thermodynamic  look  at  an  oscillating  bubble,  one 
must  examine  the  internal  pressure  of  a  bubble  as  a  function  of  time. 
Recall  that  it  was  previously  assumed  that  the  internal  pressure  is  a 
function  of  time  only.  It  has  been  predicted  that  internal  pressures  of 

hundreds  of  atmospheres  may  be  attained  during  the  violent  collapse  of  a 

bubble.  It  will  be  shown  that  under  the  proper  conditions  the  new  theory 
predicts  that  internal  pressures  on  the  order  of  hundreds  of  atmospheres 
are  attained  by  bubbles  with  large  pulsation  amplitudes.  Figures  25-27 
were  calculated  under  the  same  conditions  as  figures  22-24.  In  each 
case  the  pressure  is  normalized  with  respect  to  atmospheric  pressure. 
Note  that  figure  27  shows  an  internal  pressure  higher  than  300 
atmospheres.  These  high  pressures  last  only  a  short  time  as  do  the  high 
temperatures  discussed  previously. 

Finally,  figure  28  is  presented  with  the  radius,  internal 

pressure,  and  center  temperature  all  on  the  same  figure.  The 

equilibrium  radius  of  the  bubble  is  50  microns,  the  driving  pressure 

amplitude  0.6  bar,  and  the  driving  frequency  0.4  times  the  linear 

resonance  frequency.  From  this  figure  it  is  easy  to  see  the  phase 

relationship  between  the  radius,  pressure,  and  temperature  of  the 

bubble.  The  arrows  on  the  figure  indicate  the  positions  of  lowest 

internal  temperature  during  the  two  cycles  shown.  This  figure  supports 

our  earlier  claim,  that  the  minimum  value  of  the  gas  temperature  occurs 

when  the  radius  of  the  bubble  is  increasing  between  R  and  its  maximum 
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a  10  micron  bubble  driven  at  0.6  bar,  and  frequencies  of 
0.4  and  0.8  times  the  linear  resonance  frequency. 
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e  28.  Radius,  internal  pressure,  and  c 
function  of  time  for  a  50  micron  bubble 
frequency  of  0.4  times  the  linear  reso 


Chapter  4 

A  Systematic  approach  to  chaotic  bubble  motion 


In  recent  years,  the  numerical  investigation  of  nonlinear  dynamical 
systems  has  revealed  many  exciting  features.  Perhaps  the  most  exciting 
discovery  of  all  was  period-doubling  bifurcation  as  a  universal  route  to 
chaos  [11].  Some  other  interesting  discoveries  are  the  coexistence  of 
attractors  with  complex  basin  structures  and  the  strange  attractors 
which  exhibit  the  property  of  sensitive  dependence  on  initial 
conditions.  A  thorough  discussion  of  these  and  many  other  features  of 
nonlinear  dynamical  systems  is  found  in  the  text  by  Guckenheimer  and 
Holmes  [54]. 

In  1983  a  number  of  researchers  discovered  another  exciting 
phenomenon  associated  with  nonlinear  dynamical  systems  [55-57].  These 
researchers  have  shown  that  a  variety  of  dynamical  systems  exhibit  the 
property  of  finite  period  doubling  sequences  merging  with  inversely 
advancing  ones  to  form  a  finite  number  of  "bubbles"  on  some  cross 
section  of  the  parameter  space.  A  consequence  of  this  merging  of  period 
doubling  sequences  is  localized  regions  of  stable  orbits  and  other 
localized  regions  of  chaotic  motion  [12]. 

Probably  the  most  thoroughly  studied  dynamical  system  is  the 
Duffing  equation  [58] 

X  +  aX  +  X  +  X3  =  bcos(wt)  .  (95) 
81 


82 


One  outstanding  feature  of  this  equation  is  that  it  exhibits  all  of  the 
properties  we  have  mentioned  above.  Parlitz  and  Lauterborn  [59]  give  a 
more  global  view  of  the  Duffing  equation  by  emphasizing  the  important 
role  of  the  nonlinear  resonances,  which  cause  many  of  the  striking 
features  mentioned  above. 

One  would  like  to  have  a  global  view  of  the  new  acoustic  cavitation 
equations,  similar  to  the  one  given  for  the  Duffing  equation  in 
reference  [59].  However,  because  of  the  complexity  of  the  system,  this 
is  very  difficult.  The  many  parameters  that  can  be  varied  in  the 
acoustic  cavitation  equations  make  a  global  analysis  virtually 
impossible.  One  can,  however,  learn  a  great  deal  about  the  system  by 
using  the  techniques  described  above.  Both  the  Rayleigh-Plesset 
equation  employing  the  polytropic  approximation,  RPP,  and  the  new  exact 
formulation,  CE,  are  examined  in  this  fashon. 

One  particularly  useful  technique  is  that  of  bifurcation  analysis 
of  the  solutions.  In  general,  one  follows  the  norm  of  a  trajectory  <|> 
(  which  is  a  solution  of  the  equations  in  the  phase  space  (<t>,<t>)  )  as  it 
changes  with  respect  to  changes  in  a  control  parameter  y.  For  the 
multi-parameter  case,  one  holds  all  the  parameters  fixed  except  for  the 
one  of  current  interest  -  in  the  case  of  interest  this  is  the  driving 
pressure  amplitude.  One  then  constructs  a  bifurcation  curve  for  that 
parameter;  that  is,  one  plots  the  locus  of  points  in  the  Poincare 
section  [60]  verses  the  control  parameter. 

In  simpler  terras,  this  means  one  plots  the  value  of  the  radius  at 
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times  that  are  multiples  of  the  period  of  the  driving  pressure.  This 
technique  is  used  at  each  increment  of  the  pressure  amplitude  to 
complete  the  figure.  If  the  bubble  is  oscillating  with  a  period  equal 
to  that  of  the  driving  pressure,  one  gets  a  single  point  for  each  value 
of  the  pressure  amplitude.  If  the  period  of  the  bubble's  oscillations 
is  twice  that  of  the  driving  pressure,  two  values  of  the  radius  are 
obtained  for  each  value  of  the  pressure  amplitude,  and  so  on.  This 
type  of  graph  is  known  as  a  "Feigenbaum  tree". 

In  order  to  make  a  comparison  between  the  two  theories,  one  chooses 
the  same  set  of  initial  conditions  or  parameters.  We  choose  a  bubble 
with  an  equilibrium  radius  of  50  microns  and  a  driving  frequency  equal 
to  half  of  the  linear  resonance  frequency  of  the  bubble.  Figure  29 
shows  the  results  for  the  RPP  equation  and  figure  30  the  results  for  the 
CE  equations.  Since  the  damping  is  quite  different  for  each  model,  the 
range  of  driving  pressures  where  interesting  phenomena  occur  is 
different  as  well. 

The  most  striking  difference  between  the  two  figures,  however,  is 
that  the  RPP  equation  exhibits  a  period-doubling  bifurcation  route  to 
chaos,  while  the  CE  equations  exhibit  period  bubbling  and  an 
intermittent  transition  to  chaos.  An  interesting  feature  of  the  RPP 
Feigenbaum  tree  is  the  missing  arms  from  the  second  bifurcation  point. 
The  reason  for  the  missing  arms  is  conjectured  to  be  due  to  a  broken 
symmetry  in  the  mapping  sequence  at  the  second  bifurcation  poin'  [61]. 
From  reference  [12],  one  can  conclude  that  varying  the  driving  frequency 
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coward  a  resonance  peak  would  introduce  more  structure  in  the  figures, 
particularly  in  figure  30.  One  would  expect  to  see  period  doubling  in 
the  bubbles  of  figure  30,  and  perhaps  even  a  region  of  chaos  in  the 
bubble  itself. 

Unfortunately,  the  computer  time  required  for  a  "Feigenbaum  tree" 
using  the  CE  equations  is  over  one  week  on  a  Digital  PDP  11/73.  This 
time  constraint  has  severely  limited  the  scope  of  the  present 
investigation.  As  faster  microcomputers  become  available,  this  type  of 
analysis  is  sure  to  become  more  and  more  popular. 

Although  global  predictions  cannot  be  made  for  either  of  the  two 
theories  investigated,  we  have  learned  that  the  damping  used  in  each 
system  determines,  to  a  large  degree,  the  characteristics  of  the 
oscillations.  Thus,  one  can  compare  an  analysis  of  this  type  to 
"Feigenbaum  trees”  obtained  experimentally,  to  determine  which  model  of 
the  damping  is  more  accurate.  An  experiment  of  this  type  is  planned, 
using  the  technique  described  in  section  D  of  chapter  2. 
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Conclusions  and  topics  for  future  study 

The  purpose  of  this  study  was  to  show  that  an  exact  formulation  for 
the  internal  pressure  of  a  cavitating  bubble  leads  to  a  more  applicable 
solution  of  an  equation  modeling  acoustic  cavitation.  The  strongest 
evidence  obtained  to  support  this  claim  is  the  reasonably  close 
agreement  between  experimental  levitation  numbers  and  the  ones  obtained 
numerically  using  the  new  theory.  The  exact  formulation  shows  the 
strong  second  harmonic  resonance  peak  in  the  levitation  number  curves, 
whereas  the  polytropic  approximation  shows  very  little  resonance  effect 
at  all. 

One  can  conclude  from  this  study,  that  for  small  pulsation 
amplitudes  and  driving  frequencies  below  one  fourth  of  the  resonance 
frequency  of  a  bubble,  the  polytropic  approximation  is  sufficiently 
accurate  and  should  be  used  when  computer  time  is  a  factor.  Under  other 
conditions,  however,  the  new  formulation  should  be  used  unless  the 
pulsation  amplitude  is  large.  For  large  pulsation  amplitudes  one  cannot 
be  so  positive  about  the  new  theory.  Some  of  the  approximations  are  no 
longer  valid  when  the  pulsation  amplitude  of  the  bubble  approaches  twice 
the  equilibrium  radius.  One  must  wait  for  new  experimental  results 
before  saying  how  accurate  the  new  theory  is  in  this  region. 

It  has  been  observed  that  the  polytropic  approximation  gives  a 
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value  for  the  damping  that  is  too  large  when  the  driving  frequency  is 
near  a  harmonic  resonance.  There  are  a  number  of  ways  to  modify  the 
damping  term  in  equations  using  the  polytropic  approximation  to  help 
compensate  for  the  overdamping  at  the  resonance  peaks.  Some  of  these 
methods  may  give  results  as  accurate  as  the  exact  method.  However,  this 
is  not  very  pleasing  to  the  theoretician  interested  in  understanding  the 
physics  of  a  cavitation  bubble.  Unfortunately,  one  must  pay  a  high 
price  for  this  elegant  theory  in  terms  of  the  complexity  of  its 
numerical  solution.  The  exact  formulation  requires  about  an  order  of 
magnitude  more  computer  time  than  does  the  polytropic  formulation  with 
an  artificial  damping  term  included. 

The  artificial  damping  techniques  have  some  difficulties  of  their 
own,  however.  There  is  no  one  corrective  term  of  this  type  that  works 
well  over  a  large  range  of  equilibrium  radii,  driving  frequencies,  and 
driving  pressure  amplitudes.  This  is  a  serious  limitation  when  trying 
to  predict  the  global  behavior  of  an  acoustic  cavitation  theory. 

One  possible  solution  to  this  problem  is  to  replace  the  polytropic 
exponent  and  artificial  damping  terms  with  a  single  function  for  the 
internal  pressure.  As  stated  previously,  this  function  for  the  internal 
pressure  must  be  dependent  on  the  equilibrium  radius,  the  driving 
frequency,  and  the  driving  pressure  amplitude.  A  study  such  as  the  one 
here  may  be  of  great  help  in  devising  such  a  function.  Figures  have  been 
presented  for  the  internal  pressure  at  various  equilibrium  radii, 
driving  frequencies,  and  driving  pressure  amplitudes.  These  figures  are 
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very  useful  when  trying  to  find  the  functional  dependence  of  the 
internal  pressure  on  the  parameters  mentioned. 

In  several  of  the  figures  presented  in  the  previous  chapters,  the 
new  theory  was  pressed  to  the  limits  of  its  applicability.  One  would 
like  to  extend  the  applicabiblity  of  the  exact  formulation  into  the 
region  of  large  pulsation  amplitudes.  There  are  essentially  two 
problems  which  must  be  overcome  in  order  to  extend  the  theory  to  this 
region. 

First,  the  second  order  numerical  integrator  needs  to  be  replaced 
by  the  fourth  order  Adams-Moulton  method.  Also,  one  needs  to  find  a 
more  accurate  method  than  the  Crank-Nicolson  method  used  to  solve  for 
the  internal  temperature.  These  higher  order  methods  would  greatly 
reduce  the  numerical  error  made  during  the  collapse  phase  of  the 
bubble's  cycle  when  the  pulsation  amplitude  is  large. 

The  second  problem  encountered  is  that  of  inapplicability  of  the 
equations  for  the  internal  pressure  when  the  pulsation  amplitude  of  the 
bubble  is  too  large.  Recall  that  in  deriving  these  equations,  it  was 
assumed  in  several  places  that  the  velocity  of  the  bubble  wall  was  small 
compared  to  the  speed  of  sound  in  the  gas.  However,  this  condition  does 
not  hold  when  the  pulsation  amplitude  of  the  bubble  is  very  large. 
Hence,  one  must  retain  some  of  the  terms  discarded  in  the  earlier 
derivation  in  order  to  maintain  applicability  in  the  large  pulsation 
amplitude  region.  This  will  in  turn  make  the  numerical  integration  more 
complex  and  time  consuming  on  the  computer.  Further  investigations  are 
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needed  to  find  the  most  useful  model  which  takes  into  account  all  of  the 
difficulties  mentioned  above. 

In  previous  chapters  several  experiments  were  described  that  would 
help  to  advance  the  theory  of  acoustic  cavitation.  The  importance  of 
these  experiments  cannot  be  overemphasized.  Experimental  data  at  large 
pulsation  amplitudes  that  is  sensitive  to  the  damping  is  very  important 
to  the  theory. 

Lauterborn  has  carried  out  several  experiments  in  which  he  examines 
the  acoustic  emissions  of  a  bubble  field  and  has  been  able  to  detect  the 
onset  of  different  modes  of  oscillation  as  the  driving  pressure  is 
increased.  The  onset  of  these  different  modes  corresponds  to  a 
bifurcation  point  on  a  "Fiegenbaum  tree".  An  experiment  of  this  type 
could  be  conducted  with  a  single  bubble  in  a  levitation  cell,  such  as 
the  one  used  by  Crum.  This  would  give  an  experimental  "Feigenbaum  tree" 
which  could  be  used  to  test  the  theories.  In  chapter  4  it  was  seen  that 
the  bifurcation  points  were  very  sensitive  to  the  damping,  so  this  is  an 
appropriate  test  for  the  damping  models. 

The  laser  light  scattering  experiment  described  in  chapter  2  is 
probably  the  most  important  experiment  of  all  the  ones  mentioned.  This 
experiment  could  provide  an  exact  radius  verses  time  curve  that  would  be 
very  useful  in  developing  an  appropriate  model  of  the  internal  pressure 
or  damping. 

Until  the  successful  completion  of  some  of  the  experiments 
suggested,  or  perhaps  some  others  not  mentioned  here,  one  can  only 
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speculate  as  to  the  applicability  of  the  present  theory  to  bubbles  with 
large  pulsation  amplitudes.  One  can,  however,  say  that  this  exact 
approach  to  the  internal  pressure  of  a  cavitating  bubble  gives  more 
accurate  results  in  its  region  of  applicability  than  does  the  polytropic 
approximation.  It  is  also  pleasing  to  the  physicist  to  have  a 
theoretical  model  that  agrees  with  experiments  without  resorting  to 
artificial  viscosities  or  damping  terms. 
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APPENDIX 


A.  Computer  program  for  the  exact  formulation 

The  following  FORTRAN  program  was  used  to  make  calculations  for  the 
compressible  exact  theory,  which  was  denoted  by  CE.  This  program 


outputs 

the 

radius , 

internal  pressure, 

and  center 

temperature  as  a 

function 

of 

time. 

The  modif ications 

necessary 

to  calculate  the 

levitation  numbers,  frequency  response  curves,  temperature  profiles,  and 
Feigenbaum  tree  are  not  shown.  However,  each  of  these  modifications  is 
straightforward  to  implement  in  the  computer  program. 

The  input  parameters  necessary  to  run  the  program  are  REQ,  AMP, 
LAMBDA,  NN,  NTMSTP ,  NCYCLE ,  AND  IPRINT,  where 
REQ  is  the  equilibrium  radius  of  the  bubble  in  microns, 

AMP  is  the  driving  pressure  amplitude  in  atmospheres, 

LAMBDA  is  the  fraction  of  the  resonance  frequency  at  which  the  bubble  is 
driven, 

NN  is  the  number  of  nodes  (  finite  difference  points  )  in  the  bubble, 
NTMSTP  is  the  number  of  time  steps  in  1  period  of  the  driving  pressure, 
NCYCLE  is  the  number  of  cycles  of  the  driving  pressure  integrated  over, 
IPRINT  is  the  fraction  of  the  points  that  are  output  (  1/IPRINT  ). 
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PROGRAM  CE 


VERSION  1.0  MARCH  5,  1985 

THIS  PROGRAM  USES  A  PREDICTOR-CORRECTOR  METHOD  TO  SOLVE  THE 
NONLINEAR  SET  OF  EQUATIONS  FOR  RADIAL  BUBBLE  OSCILLATIONS. 


****************************************************************** 


DEFINITION  OF  CONSTANTS  &  VARIABLES 


C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

c 

c 

c 

c 

c 

c 

c 


GAM  =  RATIO  OF  SPECIFIC  HEATS  -  DIMENSIONLESS 
RHOL  =  DENSITY  OF  LIQUID  -  GRAMS/CC 

F  =  DRIVING  FREQUENCY  -  HERTZ 

REQ  =  EQUILIBRIUM  RADIUS  OF  BUBBLE  -  CENTIMETERS 
PINF  =  UNDISTURBED  LIQUID  PRESSURE  -  DYNES/CM/CM 
SIGMA  =  SURFACE  TENSION  -  DYNES/CM 
AMU  =  LIQUID  VISCOSITY  -  POISE 

AMP  -  DRIVING  PRESSURE  AMPLITUDE  IN  ATMOSPHERES 
EPS  =  DIMENSIONLESS  PRESSURE  AMPLITUDE  AMP/PINF 
KINF  =  THERMAL  CONDUCTIVITY  OF  AIR  @  EQUILBRIUM  COND.  - 
ERGS/ (SEC*CM*DEG. KELVIN) 

TINF  *  AMBIENT  LIQUID  TEMPERATURE  -  DEGREES  KELVIN 

PO  =  INTERNAL  PRESSURE  OF  BUBBLE  @  EQUIL.  -  DYNES/CM/CM 

W  =  RADIAL  DRIVING  FREQUENCY  -  RADIANS/SEC 

DO  =  THERMAL  DIFFUSIVITY  OF  AIR  @  EQUIL.  -  CM*CM/SEC 

NN  -  NUMBER  OF  FINITE  DIFFERENCE  POINTS  IN  BUBBLE 

NTMSTP=  NUMBER  OF  DIMENSIONLESS  TIME  STEPS  IN  ONE  PERIOD 

R  =  DIMENSIONLESS  RADIUS 

RTL  =  R  TILDA  @  N+l  TIME  STEP 

RP1  =  DIMENSIONLESS  @  N+l  TIME  STEP 

RST  =  R  STAR 

U  =  DIMENSIONLESS  VELOCITY  (BUBBLE  WALL) 

UTL  =  U  TILDA  @  N+l  TIME  STEP 

UPR  =  DU/DT  @  TIME  STEP  N 

UP1  =  DIMENSIONLESS  VELOCITY  @  N+l  TIME  STEP 

UST  =  U  STAR 

P  =  DIMENSIONLESS  INTERIOR  PRESSURE 

PTL  =  P  TILDA  @  N+l  TIME  STEP 

PPR  =  DP/DT  @  TIME  STEP  N 

PTLPR  =  DPTL/DT  @  N+l  TIME  STEP 

PST  =  P  STAR 

PSTPR  =  P  STAR  PRIME 

PP1  =  DIMENSIONLESS  INTERIOR  PRESSURE  @  N+l  TIME  STEP 
T  =  DIMENSIONLESS  TIME 

K  =  DIMENSIONAL  THERMAL  CONDUCTIVITY  OF  AIR 

VEL  =  VELOCITY  OF  SOUND  IN  LIQUID  IN  CM/SEC 
D  =  DIMENSIONLESS  THERMAL  DIFFUSIVITY  OF  AIR 
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OMO=DSQRT(OMO) 

DO  20  K-1,20 

THETA=R*DSQRT(2.*OMO/DIFF) 

EX=0 .0 

IF(  THETA . LT . 1 0 . 0 )EX=DEXP ( -THETA  ) 
ST*=2.*DSIN(  THETA) 

CT=2.*DC0S( THETA) 

APLUS=( 1 .+EX*( ST-EX) )/( 1 .+EX*(EX-CT) ) 
AMINUS=*(  1  .-EX*(ST+EX)  )/(  1  .+EX*(  EX  -CT)  ) 
T1*(THETA+THG1*AMIN(JS)*TH£  "A 
T2-THGl*(THETA*APLUS-2. ) 

DEN=T1*T1+T2*T2 
KAPPA=GAM*THETA*THETA*T1  /DEN 
OMOLD-OMO 

OMO»PRR2*(3.*KAPPA+(3.*KAPPA-1 . )*W) 
OMO=DSQRT( OMO ) 

IF( DABS( OMO/OMOLD-1 .0 ) . LT . I .D-4)C0T0  40 
20  CONTINUE 
GOTO  400 
40  CONTINUE 

AIMPHI-3 .*GAM*T2*THETA*THETA/DEN 
BETTH«.5*PIO*AIMPHI/(RHOL*OMO*R*R) 
BETAC- . 5*OMO*OMO*R/VEL 
BETAV-2 - *AMUL/ ( RH0L*R*R) 

FO=OMO/6. 283185307  179586 
THETA=R*  DS  QRT ( 2 . *  OMO / D I FF ) 
BETA=BETTH+BETAV+BETAC 
400  CONTINUE 
RETURN 
END 
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DDR2=H*DDER 

R=RI+.5*DR2 

DR=DRI+.5*DDR2 

DDER=DER(R,DR,TT) 

DR3=H*DR 

DDR3=H*DDER 

R=RI+DR3 

DR=DRI+DDR3 

TT=T+H 

DDER=DER( R , DR , TT ) 

DR4=H*DR 

DDR4=H*DDER 

R0=RI+( DR1+DR4+2 . * ( DR2+DR3 ) ) / 6 . 
DR0=DRI+(DDRl+DDR4+2 .*( DDR2+DDR3 ) )/6 . 

RETURN 

END 

DOUBLE  PRECISION  FUNCTION  DER( R, DR,T) 

IMPLICIT  REAL* 8 ( A-H , O-Z ) 

REAL* 8  KAPPA, LAMBDA 
COMMON/ SET 1 /AMP , KAPPA , LAMBDA 
COMMON/ COST /PAR , AMTH , AMEFF , P I , P IPOL 
COMMON/ Q/W,PI0 
COMMON/SET4 /THREEK , OMW 

PIPOL  *  P  INTERNAL  POLYTROPIC 

PI  *  PIPOL  -  THERMAL  DAMPING  CONTRIBUTION 

PIPOL=l./R**THREEK 

PI=PIPOL-AMTH*DR/R 

DER»PI-W/R-AMEFF*DR/R-OMW*( 1 .-AMP*DSIN(T) ) 

DER»(PAR*DER-1.5*DR*DR)/R 

RETURN 

END 

SUBROUTINE  RESNCE(R) 

IMPLICIT  REAL*8(A-H,0-Z) 

REAL* 8  KAPPA 

COMMON/ SET 1 /AMP , KAPPA , LAMBDA 

COMMON/PARAMS/GAM , DIFF , PINF , RHOL , SIGMA , VEL , AMUL 
COMMON/RES/ FO , AIMPHI , BETTH , BETAV , BETAC , BETA 
COMMON /Q/W,P 10 
PI0-PINF+2.*SIGMA/R 
THG1*3.*( GAM-1.) 

PRR2«PINF/(RH0L*R*R) 

W«2.*SIGMA/(R*PINF) 

KAPPA=1 .2 

OMO=PRR2*(3.*KAPPA+(3.*KAPPA-l . )*W) 


1 

1 


i 


H2=H/2. 

EPS=l.D-7 

EPSB=l.D-8 


DO  40  K- 1,32500 

CALL  RUNGE(RO,DRO,RA,DRA,T,H) 

20  CONTINUE 

CALL  RUNGE ( RO , DRO , RB , DRB , T , H2 ) 

T=T+H2 

CALL  RUNGE(RB,DRB,RC,DRC,T,H2 ) 

T=*T+H2 

IF(DABS(DRC).LT.EPSB)GOTO  30 
ERR«DABS( DRA/DRC-1 . ) 

IF( ERR.LT.EPS)  GO  TO  30 

RA=RB 

DRA=DRB 

T=T-H 

H=»H2 

H2-H/2. 

GO  TO  20 
30  CONTINUE 
R0=RC 
DR0=DRC 

DDER=DER(R0,DR0,T) 

35  PRES=DSIN(T) 

WRITE( 1 ,38 )R0  ,T 
WRITE (2 ,38)PI,T 
IF(T.GE.TMAX)GOTO  200 
38  FORMAT(F12.6) 

36  IF(ERR.GT.EPSB)GO  TO  40 
H2*H 
H-2.*H2 
40  CONTINUE 
200  CONTINUE 
RETURN 
END 
C 

SUBROUTINE  RUNGE(RI , DRI , RO , DRO ,T , H) 
IMPLICIT  REAL* 8 ( A-H , 0-Z ) 
DDER-DER(RI.DRI.T) 

DR1=*H*DRI 
DDR1-H*DDER 
R=»RI+.5*DR1 
DR=DRI+ • 5*DDR1 
TT-T+H/2. 

DDER*DER(R,DR,TT) 

DR2=H*DR 
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COMMON/PARAMS/GAM.DIFF.PINF, RHOL, SIGMA, VEL.AMUL 

BETTHR=BETTH 

AIMPHR=AIMPHI 

BETACR=BETAC 

R0=1 .0 

DR0=*0 .0 

H-6. 283185307 179586D-2 
OMW=l .-W 

OMO«=6. 283185307 179586*F0 
OM=  LAMB  DA*  OMO 

PAR=P 1 0 / ( RHOL* ( 0M*R ) * ( 0M*  R ) ) 

AMVIS-4 . *AMUL*0M/PI0 
THETA«R*DSQRT(2.*0M/DIFF) 

EX=*0 .0 

IF(THETA.LT.10.0)EX=DEXP(-THETA) 

ST-2.*DSIN(THETA) 

CT-2.*DC0S(THETA) 

APLUS*( 1 .+EX*(ST-EX) )/( 1 .+EX*(EX-CT) ) 

AMINUS= ( 1 .-EX* ( ST+EX ) ) / ( 1 . +EX* ( EX-CT ) ) 

T 1 - ( THETA+3 . * ( GAM-1 ) * AMINU  S ) * THETA 
T2-3 . * ( GAM-1 . ) * ( THETA* APLUS-2 . ) 

DEN«T1*T1+T2*T2 

AIMPHN-3 .*GAM*T2*THETA*THETA/DEN 

KAPPA=GAM*THETA*THETA*T 1 /DEN 

THREEK- 3.0* KAPPA 

BETACN- . 5*0M*0M*R/VEL 

BETTHN- . 5*PI0*AIMPHN/ ( RHOL*R*R*OM ) 

BETTH-BETTHR 
AIMPHI-AIMPHR 
BETAC-BETACR 
IF(OPTTH.EQ. 1 .0 )G0T0  10 

NOW  COMPUTE  THE  THERMAL  DAMPING  AT  OM  RATHER  THAN  OMO 

AIMPHI-AIMPHN 

BETTH-BETTHN 

BETAC-BETACN 

CONTINUE 

AMAC- 2 . *RHOL*OM* R*  R*  BETAC /P 10 
AMTH-AIMPHI 

NEXT  STATEMENT  IS  FOR  NO  THERMAL  DAMPING 
AMTH-0.0 


AMEFF-AMVIS+AMAC 

TMAX=6. 283185307 179586*NNNN 

T-0.0D0 


o  o 


PROGRAM  RPP 

IMPLICIT  REAL*8 ( A-H , O-Z ) 

REAL* 8  KAPPA, LAMBDA, RDSK 
COMMON/COST/PAR, AMTH,AMEFF, PI, PIPOL 
COMMON/ SET 1 /AMP  .KAPPA , LAMBDA 
COMMON/ SET2/NNNN 
COMMON/ SET3/WRK, NTM1 ,DT 

COMMON/PARAMS/GAM , DIFF , PINF , RHOL , S IGMA , VEL , AMUL 

OPTTH-O.O  GIVES  NONRESONANCE  VALUES: 1.0  GIVES  RESONANCE 
0PTTH=0 .0 
C 

GAM= 1.4 
DIFF-0 .2 
PINF=1 .013D6 
RHOL- 1.0 
AMUL-1.0D-2 
SIGMA-72.8 
VEL-1.481D5 
C 

WRITE(5 ,55) 

55  FORMAT ( '  INPUT  RO,  LAMBDA,  PA,  NO.  OF  CYCLES') 

READ ( 5 , * ) RDIM , LAMBDA , PA , NNNN 

C 

R-RDIMM.OD— 4 
AMP- PA 
PA-PA* 1 .0D6 
C 

CALL  RESNCE(R) 

C 

WRITE (5, 56) 

56  FORMAT ( 1H  ) 

C 

CALL  WORK ( R , OPTTH ) 

C 

STOP 

END 

C 

SUBROUTINE  WORK (R, OPTTH) 

IMPLICIT  REAL*8( A-H.O-Z) 

REAL*8  KAPPA, LAMBDA 

COMMON/ SET5 / DELAY , TR1 

COMMON/SET2/NNNN 

COMMON/ SET 1 / AMP , KAPPA , LAMBDA 

COMMON/ SET4 /THREEK , OMW 

COMMON/Q/W.PIO 

COMMON/ RES / FO , AIMPHI , BETTH , BETAV , BETAC ,BETA 
COMMON/COST/PAR, AMTH.AMEFF, PI, PIPOL 
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B.  Computer  program  for  the  Rayleigh-Plesset  polytropic  equation 

The  following  FORTRAN  program  was  used  to  make  calculations  for  the 
Rayleigh-Plesset  equation  using  the  polytropic  approximation.  The 
integration  is  carried  out  using  a  fourth  order  Runge-Kutta  algorithm. 
The  output  from  this  program  is  radius  and  internal  pressure  as  a 
function  of  time. 

The  input  parameters  necessary  to  run  the  program  are  RDIM, 

LAMBDA,  PA,  NNNN,  where, 

RDIM  is  the  equilibrium  radius  of  the  bubble  in  microns. 

Lambda  is  the  fraction  of  the  resonance  frequency  at  which  the  bubble  is 
driven, 

PA  is  the  driving  pressure  amplitude  in  atmospheres,  and 

NNNN  is  the  number  of  cycles  of  the  driving  pressure  amplitude 


integrated  over 


o  o  o 


105 


RETURN 

END 

TRIDIAGONAL  SYSTEM  SOLVER 
SUBROUTINE  TRIDG(VEC) 

DIMENSION  Dll (5 1 ) ,DI2(51 ),DI3(51), RHS(51 ) , VEC(52 ) 
COMMON/SET2/DI1 ,DI2 ,DI3 , RHS , NM1 ,N 
DO  10  I-l.NMl 

DI2(I+1)-DI2(I+1)“DI3(I)*DI1(I)/DI2(I) 

10  RHS(I+1)-RHS(I+1)-8HS(I)*DI1(I)/DI2(I) 

VEC(N)=RHS(N)/DI2(N) 

DO  20  M=NM1 , 1 ,-l 

20  VEC(M)=(RHS(M)-DI3(M)*VEC(M+1))/DI2(M) 

RETURN 

END 


o  n  n  no 
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WRITE( 1 ,96)R,T 
WRITE(2,96)P,T 
WRITE ( 3,96)TG(1),T 
96  FORMAT(F12.5) 

1000  CONTINUE 
1100  CONTINUE 


STOP 
END 

SUBROUTINE  TO  CALCULATE  RESONANCE  FREQUENCY  FO 

SUBROUTINE  RESNCE(R.FO) 

REAL  KAPPA, LAMBDA 

COMMON/ SET 1 /GAM , P I NF , RHOL , S IGMA , VEL , AMUL , KINF , T I NF 
P IO-P INF+2 . *  S IGMA/R 
THG1-3.*(GAM-1.) 

DIFF-(GAM-1.0)*KINF*TINF/(GAM*PI0) 

DIFF-.2 

PRR2-PINF/ ( RHOL*R*R ) 

W-2. *S IGMA/ (R*P INF) 

KAPPA- 1.2 

OMO«PRR2*(3 .*KAPPA+(3 .*KAPPA-1 . )*W) 

OMO-SQRT(OMO) 

DO  20  K-1,20 

THETA=R*SQRT(2.*OMO/DIFF) 

EX-0.0 

IF( THETA.LT. 10 .0 )EX-EXP( -THETA) 

ST-2 . *SIN(THETA) 

CT-2.*COS(THETA) 

APLUS-(1 .+EX*(ST-EX))/( 1 .+EX*(EX-CT) ) 

AMINUS- ( 1 .-EX* ( ST+EX ) ) / ( 1 . +EX* ( EX-CT ) ) 

T 1 - ( THETA+THG 1 *AMI NUS ) *  THET A 
T2-THGl*(THETA*APLUS-2 . ) 

DEN-T1*T1+T2*T2 
KAPPA-GAM*THETA*THETA*T 1 /DEN 
OMOLD-OMO 

OMO-PRR2*(3.*KAPPA+(3.*KAPPA-l.)*W) 

OMO-SQRT ( OMO ) 

IF(ABS(OMO/OMOLD-1.0).LT.l .E-4)G0T0  40 
20  CONTINUE 

WRITE (5 ,30) 

30  FORMATC  NO  CONVERGENCE  IN  SUBROUTINE  RESNCE ' ) 
GOTO  400 
40  CONTINUE 

F0-0M0/6. 2831853 
400  CONTINUE 


ooono  ooo 
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C 

RST2=RST*RST 

DUM1-3.0*DTDY2*CHI*DST( 1)/RST2 

DI2(1)-1.0+DUM1 

DX3( 1 )=-DUMl 

RHS( 1 )-TAU( 1 )+.5*DT*(DTAUDT( 1 )+DST( 1 )*PSTPR) 

DO  80  1-2, NN 

TERM-(DTDY2*CHI/(2.0*RST2))*(GM1G*(.5*(TAUST(I+1)-TAUST(  1-1)) 
1+Y( I )*TAUST( NN) ) /PST ) 

Dll ( 1-1 )— TERM-( DTDY2*CHI/( 2 .0*RST2 ) )*DST( I )*YM( I ) 

DI2( I)=l .0+DTDY2*CHI*DST( I)/RST2 
DI3(I)-TERM-DST(I)*YP(I)*(DTDY2*CHI/(2.0*RST2)) 

RHS( I)-TAU( I )+ . 5*DT*(DTAUDT( I )+DST( I )*PSTPR ) 

80  CONTINUE 
C 

CALL  TRIDG(TAUPl) 

OTHER  CORRECTOR  EQUATIONS 
MACH-UST/C 

DRIVER=PFP0*( 1 .0-EPS*SIN(T+DT+RST/DUME) ) 

DUM1-1 .0+(DT/(2 .0*RST*( 1 .0-MACH) ) )*( 1 .5*UST*( 1 .0-MACH/3 .0)+CAPP* 
1(1 .0+MACH)*CAPM/RST) 

DUM2-U+.5*DT*(UPR+CAPP*((PST-DRIVER-CAPW/RST)*( 1 .0+MACH)+ 

1 ( 1 .0+0 .0 )*RST*PSTPR/C ) /(RST*( 1 .0-MACH) ) ) 

UP1-DUM2/DUM1 

RP1-R+.5*DT*(U+UP1) 

PP1-P+ .5*DT*(PPR-3 .0*GM1*CHI*TAUP1 ( NN)/( DY*RP1*RP1 ) ) 

PP1-PP1/ (1.0+1 .5*GAM*DT*UP1 /RP1 ) 

DO  85  1-1 , NN 

TGP 1 ( I )- ( SQRT( 1 .0+GM1 G2A*TAUP1 ( I ) )-BETA ) /ALPHA 
DP1( I)-(ALPHA*TGP1( I)+BETA)*TGP1 ( I)/PP1 
85  CONTINUE 

****************************************************************** 


UPDATE  AND  OUTPUT  AREA 

T-TCC+FLOAT (J)*6. 28318531/FLO AT (NTMSTP) 
R-RP1 
P-PP1 
U-UP1 

DO  90  I- 1 , NN 
TG( I)-TGP1( I) 

TAU(I)-TAUPKI) 

D(I)-DP1(I) 

90  CONTINUE 

IF(M0D( J, IPRINT) .NE.O)GOTO  1000 


n  n  o  o 


► 
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c 


« 


r 


( 


RTL=R+DT*U 

RTL2=RTL*RTL 

DRIVER*PFP0*( 1 .0-EPS*SIN(T+R/DUME) ) 

UPR— 1 . 5*U*U* ( 1 . -MACH/3 . )+CAPP* ( ( P-DRIVER-( CAPW+CAPM*U ) /R) 

1*( 1 .0+MACH)+( 1 .0+0 .0 )*R*PPR/C )/( R*( 1 .0-MACH) ) 

UTL=U+DT*UPR 

PTLPR*3.0*(GM1*CHI*(-TAUTL(NN))/RDY-GAM*PTL*UTL)/RTL 

c 

c  START  OF  FIRST  CORRECTOR  SECTION 

C 

C  TRI-DIAGONAL  SOLUTION  FOR  TAUST  *  TAU  STAR  @  N+l  STEP 

C 

DUM1*3 .0*DTDY2*CHI*DTL( 1 )/RTL2 
DI2( 1 )=1 .0+DUM1 
DI3( 1 )*-DUMl 

RHS ( 1 ) =TAU ( 1 )+DT* ( DTAUDT ( 1 )+ DTL( 1 ) * PTLPR ) / 2 . 0 
DO  70  1*2, NN 

TERM-(DTDY2*CHI/(2.0*RTL2))*(GM1G*( .5*(TAUTL(I+1 )-TAUTL( 1-1 ) ) 
1+Y( I)*TAUTL( NN) )/PTL) 

Dll  ( 1-1 )— TERM-(  DTDY2*CHI/  (  2  .0*RTL2  )  )*DTL(  I  )*YM(  I ) 

DI2( I)*l .0+DTDY2*CHI*DTL( I)/RTL2 

DI3 ( I )=TERM-DTL( I)*YP( I )*( DTDY2 *CHI/ ( 2 . 0*RTL2 ) ) 

RHS( I )*TAU( I )+ . 5*DT*( DTAUDTC I )+DTL( I )*PTLPR) 

70  CONTINUE 

C 

CALL  TRIDG( TAUST) 

C 

C  OTHER  CORRECTOR  EQUATIONS 

C 

MACH=UTL/C 

DUM1*1 .0+(DT/(2 .0*RTL*( 1 .0“ MACH) ))*( 1 .5*UTL*( 1 .O-MACH/3 .0)+CAPP* 
1(1.0+MACH)*CAPM/RTL) 

DRIVER*PFP0*( 1 .0-EPS*SIN(T+DT+RTL/DUME) ) 

DUM2*U+ . 5*DT* ( UP R+CAP  P* ( ( PTL-DR IVE  R-CAPW/RTL ) * ( 1 . 0+MACH ) + 

1 ( 1 . 0+0 . 0 ) *RTL* PTLPR/ C ) / ( R* ( 1 . 0 -MACH ) ) ) 

UST*DUM2/DUM1 
RST=R+ . 5  *  DT* ( U+UST ) 

PST*P+ . 5*DT* (PPR-3 . 0*GM1 *CHI*TAUST( NN) / ( DY*RST*RST ) ) 

PST*PST/( 1 .0+1 ,5*GAM*DT*UST/RST) 

PSTPR*3 .0*(GM1*CHI* (-TAUST ( NN) )/( RST*DY)-GAM*PST*UST)/RST 
DO  75  I-l.NN 

TGST( I)-(SQRT( 1 .0+GMlG2A*TAUST( I ) )-BETA) /ALPHA 
DST(I)-(ALPHA*TGST(I)+BETA)*TGST(I)/PST 

75  CONTINUE 

SECOND  CORRECTOR  EQUATIONS  START  HERE 
TRI-DIAGONAL  SOLUTION  FOR  TAUP1  *  NEW  VALUE  OF  TAU  @  N+l  STEP 


4 


non  o  n  n  onnoo  ooo 
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Y(II)=Y(II-1)+DY 

YP(II)=1.0+DY/Y(II) 

YM(II)=1.0-DY/Y(II) 

TG(II)=1.0 
TAU(II)=0.0 
50  CONTINUE 
C 

£*★*★** ******* ******** **************  ********** ************* ***  ****** 


USEFUL  EXPREESSIONS  TO  SAVE  COMPUTATIONAL  TIME 

DUME=W*C 
DY2=DY*DY 
DY22=DY2*2 .0 
DTDY2=DT/DY2 


WORK  AREA 


DO  1100  11=1 ,NCYCLE 
TCC=(II-1)*6. 283185307 
TYPE  *,II 

DO  1000  J= 1 , NTMSTP 

START  OF  PREDICTOR  ROUTINE 


R2=R*R 

RDY=R*DY 

R2DY2=R2*DY2 

PPR=-3 . 0*( GM1 *CHI*TAU(NN ) /RDY+GAM*P*U ) /R 
PTL=P+DT*PPR 

DTAUDTC 1 )»D( 1 )*(6.0*CHI*(TAU(2 )~TAU( 1 ) )/R2DY2+PPR) 

TAUTL( 1 )=TAU( 1 )+DT*DTAUDT( 1 ) 

TGTL( 1 )=(SQRT( 1 .0+GM1 G2A*TAUTL( 1 ) )-BETA) /ALPHA 
DTL( 1 )=( ALPHA*TGTL( 1 )+BETA)*TGTL( 1 ) /PTL 
DO  60  1=2, NN 

DTAUDT( I)=D( I)*( (CHI/R2DY2 )*( YP( I )*TAU( 1+1 )-2 .0*TAU( I)+ 

1YM( I)*TAU( 1-1 ) )+PPR)-(GMlG*CHl/( R2*P) )*( (TAU( 1+1 )-TAU( I-l))/2.0 
2+Y( I)*TAU( NN) )*( TAU( 1+1 )-TAU( 1-1 ) ) /DY22 
TAUTL( I )»TAU( I )+DT*DTAUDT( I ) 

TGTL( I)=(SQRT( 1 .0+GMlG2A*TAUTL( I) )-BETA) /ALPHA 
DTL( I )=( ALPHA*TGTL( I)+BETA)*TGTL( I ) /PTL 
CONTINUE 

REST  OF  PREDICTOR  CALCULATIONS 


MACH=U/C 


ooooo  noooo 


CALCULATION  OF  THE  RESONANCE  FREQUENCY  FO 
CALL  RESNCE( REQ , FO ) 

DERIVED  CONSTANTS 


F-FO*LAMBDA 

EPS-AMP/PINF 

DT-6. 2831 853 /FLOAT(NTMSTP) 

DELAY-6 .2831853*( NCYCLE-2 ) 

DY=1.0/FL0AT(NN) 

NP1-NN+1 
NM1-NN— 1 

PO-PINF+2 .0*SIGMA/REQ 
PFPO-PINF/PO 
W=6.2831853*F 
C=VEL/(W*REQ) 

GM1 -GAM-1.0 
GM1G-GM1/GAM 
D0=GM1G*KINF*TINF/P0 
CHI-DO /(W* REQ* REQ) 

CAPP=P0/( RHOL*W*W*REQ*REQ) 

CAPW-2 . 0*SIGMA/ ( REQ*P0 ) 

CAPM-4 .0*AMU*W/P0 
ALPHA=M*T INF/K INF 
BETA-B/KINF 
GM1G2A-2 ,0*GM1G* ALPHA 

************ ****** ************ ************** ******************** 


INITIAL  CONDITIONS 


R-1.0 

P-1.0 

T-0.0 

U-0.0 

D( 1 )-l  .0 

TG( I )-1 .0 

TG(NP1 )-l  .0 

TAU( 1 )-0 .0 

TAU(NPl)-0.0 

DTAUDT(NP1 )-0.0 

Y( 1 )-0.0 

Y(NP1 )-l  .0 

C  NOTE-  YP  &  YM  @  1  ARE  NOT  REQUIRED  IN  COMPUTATIONS. 
DO  50  II-2.NN 
D(II)-1.0 


noonn  non 
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V 


c 


c 

DST 

32 

D  STAR 

c 

DTL 

SC 

D  TILDA 

c 

DPI 

32 

D  ARRAY  AT  NEW 

TIME  STEP 

c 

TG 

32 

DIMENSIONLESS 

INTERIOR 

GAS 

c 

TGTL 

32 

TG  TILDA 

c 

TGST 

= 

TG  STAR 

c 

TGP1 

* 

TG  @  N+l  TIME 

STEP 

c 

TAU 

S 

DIMENSIONLESS 

INTERIOR 

GAS 

c 

TAUTL 

32 

TAU  TILDA 

c 

TAUST 

S 

TAU  STAR 

c 


TEMPERATURE  ARRAY 


TEMPERATURE  ARRAY 


£**** ** **********************************************************  ** 


c 

REAL  KINF,M,K, LAMBDA, MACH 
DIMENSION  TGTL( 51),TGST(51) 

DIMENSION  Dll (51 ) ,DI2( 51 ) ,DI3( 5 1 ) ,RHS(51 ) ,DP1 (51 ) 

DIMENSION  D(51 ) ,DST( 51 ) ,DTL( 51 ) , Y(51 ) ,YM( 51 ) ,YP( 51 ) ,DTAUDT( 51 ) 
DIMENS ION  TG( 5 1 )  , TGP 1 ( 5 1 ) , TAU( 5 1 ) ,TAUTL( 5 1 ) , TAUST( 5 1 ) , TAUP 1(51) 
COMMON/ SET1 /GAM , PINF , RHOL , S IGMA , VEL , AMU , KINF , TINF 
COMMON/SET2 /DI 1 , DI2 , DI3 , RHS , NM1 , NN 
C 

Qiticjt'kirkit'kjcjrk  ******************************************  ************* 


INPUT  FROM  TERMINAL 


WRITE(5 , 10) 

10  FORMATS  INPUT  RADIUS  IN  MICRONS,  PA,  AND  LAMBDA’) 

READ( 5 ,*) REQ , AMP , LAMBDA 
WRITE(5 ,20) 

20  FORMAT( ’  INPUT  NN,  NTMSTP,  NCYCLE ,  AND  IPRINT') 

READ( 5 ,* )NN, NTMSTP , NCYCLE , IPRINT 
REQ=REQ*1 .OE-4 
AMP=*AMP*1 .013E6 

***************************************************************** 


CONSTANTS 


M=5 .528 
B®1 165 .0 
GAM= 1.4 

<■  RHOL=0 .998 

PINF-1.013E6 
SIGMA=72 .8 
AMU-0.01 
TINF=293 . 15 
KINF=TINF*M+B 

1  VEL-1.481E5 


